{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import linalg as LA\n",
    "\n",
    "from qiskit import(\n",
    "    QuantumCircuit,\n",
    "    ClassicalRegister, \n",
    "    QuantumRegister,\n",
    "    transpile)\n",
    "\n",
    "from qiskit.circuit.library import UnitaryGate\n",
    "from qiskit.circuit.library import efficient_su2\n",
    "from qiskit.quantum_info import Statevector\n",
    "\n",
    "from qiskit.providers.basic_provider import BasicProvider\n",
    "# from qiskit_ionq import IonQProvider, GPIGate, GPI2Gate, ZZGate\n",
    "\n",
    "# backend = BasicProvider().get_backend('basic_simulator')\n",
    "# provider = IonQProvider()\n",
    "\n",
    "from qiskit_aer import AerSimulator\n",
    "from qiskit_ibm_runtime.fake_provider import FakeWashingtonV2\n",
    "backend = AerSimulator.from_backend(FakeWashingtonV2())\n",
    "\n",
    "# backend_native = provider.get_backend(\"simulator\", gateset='native')\n",
    "\n",
    "from givens_angles import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([[0, 1], [1, 0]], dtype='complex')\n",
    "y = np.array([[0, -1j], [1j, 0]], dtype='complex')\n",
    "Id = np.identity(2)\n",
    "op = np.kron(Id, x) @ np.kron(y, Id) - np.kron(Id, y) @ np.kron(x, Id)\n",
    "\n",
    "def givens_circ(theta):\n",
    "    mat = LA.expm(1j/2 * theta * op)\n",
    "    qreg = QuantumRegister(2)\n",
    "    circ = QuantumCircuit(2)\n",
    "    circ.append(UnitaryGate(mat), [0, 1])\n",
    "    return circ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N_s = 8\n",
    "N_t = 1\n",
    "J = 0.4\n",
    "h = 1\n",
    "g = 0.3\n",
    "dt = 0.6\n",
    "N_t * dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_op(N_s, site, op):\n",
    "    Op = np.identity(1)\n",
    "    for s in range(N_s):\n",
    "        if s == site:\n",
    "            Op = np.kron(Op, op)\n",
    "        else:\n",
    "            Op = np.kron(Op, np.identity(2))\n",
    "    return Op\n",
    "\n",
    "x = np.array([[0, 1], [1, 0]])\n",
    "y = np.array([[0, -1j], [1j, 0]])\n",
    "z = np.array([[1, 0], [0, -1]])\n",
    "n = np.array([[0, 0], [0, 1]])\n",
    "sm = (x - 1j * y) / 2\n",
    "\n",
    "X = []\n",
    "Z = []\n",
    "N = []\n",
    "Sm = []\n",
    "\n",
    "for i in range(N_s):\n",
    "    X.append(build_op(N_s, i, x))\n",
    "    Z.append(build_op(N_s, i, z))\n",
    "    N.append(build_op(N_s, i, n))\n",
    "    Sm.append(build_op(N_s, i, sm))\n",
    "\n",
    "H_XX = np.zeros((2**N_s, 2**N_s), dtype=complex)\n",
    "H_Z = np.zeros((2**N_s, 2**N_s), dtype=complex)\n",
    "H_ZZ = np.zeros((2**N_s, 2**N_s), dtype=complex)\n",
    "\n",
    "for j in range(N_s):\n",
    "    H_XX -= J * X[j] @ X[(j + 1) % N_s]\n",
    "    H_Z -= h * Z[j]\n",
    "    H_ZZ -= g * Z[j] @ Z[(j + 1) % N_s]\n",
    "\n",
    "TEO = LA.expm(-1j * dt * H_Z) @ LA.expm(-1j * dt * H_XX) @ LA.expm(-1j * dt * H_ZZ)\n",
    "\n",
    "H = H_XX + H_Z + H_ZZ\n",
    "Evals, Evecs = np.linalg.eigh(H)\n",
    "ground_state = Evecs[:,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.00000000e+00-4.20429890e-17j  5.55111512e-17-2.44515335e-17j\n",
      " -5.55111512e-17-6.75547015e-17j -1.38777878e-17+2.48615523e-17j]\n",
      "[ 1.00000000e+00+5.52718700e-17j  5.55111512e-17+3.10407827e-17j\n",
      " -1.11022302e-16-2.82678355e-17j -2.77555756e-17+2.25654659e-17j]\n"
     ]
    }
   ],
   "source": [
    "n_k = 3\n",
    "sigma = 3/2\n",
    "# Finds Givens Angles for givens N_s, n_s (if localized),\n",
    "# mean position x, mean momentum number (integer), and width.\n",
    "# From givens_angles.py\n",
    "r = [-1, 0, 1, 2]\n",
    "betas_r, angles_r = giv_angles(N_s, 4, 1, n_k, sigma, r=r)\n",
    "betas_l, angles_l = giv_angles(N_s, 4, 5, -n_k, sigma, r=r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "vqe_angles = np.load('vqe-thetas-N8-J0.4-h1-g0.3.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_circuit(t_step):\n",
    "    qreg = QuantumRegister(N_s+4)\n",
    "    creg = ClassicalRegister(N_s+4)\n",
    "    circ = QuantumCircuit(qreg, creg)\n",
    "\n",
    "    # State-Prep\n",
    "    ## Ground state preparation\n",
    "    # circ.prepare_state(ground_state, range(N_s)) # Exact\n",
    "    ## VQE\n",
    "    for i in range(N_s):\n",
    "        circ.ry(vqe_angles[i], i)\n",
    "        circ.rz(vqe_angles[N_s + i], i)\n",
    "\n",
    "    for i in range(2, len(vqe_angles)//N_s-1,2):    \n",
    "        for j in range(N_s-1):\n",
    "            k = N_s - 1 - j\n",
    "            circ.cx(k-1, k)\n",
    "        for j in range(N_s):\n",
    "            circ.ry(vqe_angles[i * N_s + j], j)\n",
    "            circ.rz(vqe_angles[(i+1) * N_s + j], j)\n",
    "\n",
    "    circ.x(N_s) # Prepare control\n",
    "    circ.x(N_s+2) # Prepare control\n",
    "\n",
    "    ## Right moving wave packet\n",
    "    ### V^dag\n",
    "    for i in range(N_s//2):\n",
    "        circ.rz(betas_r[i], i)\n",
    "    for i in range(N_s//2-1):\n",
    "        j = N_s//2 - 1 - i\n",
    "        circ.append(givens_circ(angles_r[i]), [j-1, j])\n",
    "    ###\n",
    "    circ.ccx(N_s, 0, N_s+1) # Remove |1> part of system qubit\n",
    "    circ.x(0) # Excite system qubit\n",
    "    ### V\n",
    "    for i in range(N_s//2-1):\n",
    "        j = N_s//2 - 2 - i\n",
    "        circ.append(givens_circ(-angles_r[j]), [i, i+1])\n",
    "    for i in range(N_s//2):\n",
    "        circ.rz(-betas_r[i], i)\n",
    "    \n",
    "    # # Left Moving Wave Packet\n",
    "    # # V^dag\n",
    "    for i in range(N_s//2):\n",
    "        circ.rz(betas_l[i], i+N_s//2)\n",
    "    for i in range(N_s//2-1):\n",
    "        j = N_s - 1 - i\n",
    "        circ.append(givens_circ(angles_l[i]), [j-1, j])\n",
    "    ###\n",
    "    circ.ccx(N_s+2, N_s//2, N_s+3) # Remove |1> part of system qubit\n",
    "    for i in range(N_s//2): # Obey anticommutation relations\n",
    "        circ.z(i)\n",
    "    circ.x(N_s//2) # Excite system qubit\n",
    "    ### V\n",
    "    for i in range(N_s//2-1):\n",
    "        j = N_s//2 - 2 - i\n",
    "        circ.append(givens_circ(-angles_l[j]), [i+N_s//2, i+N_s//2+1])\n",
    "    for i in range(N_s//2):\n",
    "        circ.rz(-betas_l[i], i+N_s//2)\n",
    "    \n",
    "    # Trotter Evolution\n",
    "    for i in range(t_step):\n",
    "        # Transverse Field\n",
    "        for j in range(N_s):\n",
    "            circ.rz(-2 * h * dt, j)\n",
    "        # Nearest-neighbor\n",
    "        for j in range(N_s):\n",
    "            circ.rxx(-2 * J * dt, j, (j + 1) % N_s)\n",
    "        # Longitudinal Field\n",
    "        for j in range(N_s):\n",
    "            circ.rzz(-2 * g * dt, j, (j + 1) % N_s)\n",
    "    \n",
    "    # Measure\n",
    "    circ.measure(qreg, creg[::-1])\n",
    "\n",
    "    # new_circuit = transpile(circ, backend, \n",
    "    #                     basis_gates=['cx', 'id', 'rz', 'x', 'sx'], \n",
    "    #                     optimization_level=3)\n",
    "\n",
    "    return circ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "basis_gates = ['cz', 'rx', 'rz', 'rzz', 'sx', 'x']\n",
    "\n",
    "circuits = []\n",
    "gates = []\n",
    "depths = []\n",
    "for i in range(N_t):\n",
    "    circuit = make_circuit(i)\n",
    "    new_circuit = transpile(circuit, backend=backend, optimization_level=3)\n",
    "    circuits.append(new_circuit)\n",
    "    gates.append(dict(new_circuit.count_ops()))\n",
    "    depths.append(new_circuit.depth())\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[117]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "depths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[{'rz': 151, 'sx': 110, 'cx': 67, 'measure': 12, 'x': 3}]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_shots = 100000\n",
    "job = backend.run(circuits, shots=num_shots)\n",
    "counts_list = job.result().get_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [],
   "source": [
    "# filename = 'qsim-free-Ns' + str(N_s)\n",
    "# # filename = 'ionq-run-Feb10-Ns' + str(N_s)\n",
    "# filename += '-lam' + str(lam)\n",
    "# filename += '-eps' + str(hL)\n",
    "# filename += '-dt' + str(dt)\n",
    "# filename += '-shots' + str(num_shots)\n",
    "# filename += '.json'\n",
    "\n",
    "# import json\n",
    "# with open(filename, 'w') as fout:\n",
    "#     json.dump(counts_list, fout)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_list=[counts_list]\n",
    "qc_occ_nums = []\n",
    "for i in range(N_t):\n",
    "    expecz = np.zeros(N_s)\n",
    "    num_anc_0 = 0\n",
    "    for key in counts_list[i]:\n",
    "        if key[N_s+1] == '0' and key[N_s+3] == '0':\n",
    "            num_anc_0 += counts_list[i][key]\n",
    "            for j in range(N_s):\n",
    "                if key[j] == '0':\n",
    "                    expecz[j] += counts_list[i][key]\n",
    "                else:\n",
    "                    expecz[j] -= counts_list[i][key]\n",
    "    qc_occ_nums.append((1 - expecz/num_anc_0) / 2)\n",
    "qc_occ_nums = np.array(qc_occ_nums)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.bar(range(N_s), qc_occ_nums[0])\n",
    "# plt.xlabel(r'site $j$')\n",
    "# plt.ylabel(r'$\\langle N_j \\rangle$')\n",
    "# plt.title(\"(a) Initial State\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# x = np.linspace(0, N_s-1, N_s)\n",
    "# y = np.linspace(0, (N_t - 1) * dt, N_t)\n",
    "# X, Y = np.meshgrid(x, y)\n",
    "# fig, ax = plt.subplots()\n",
    "\n",
    "# plt.pcolormesh(X, Y, qc_occ_nums, cmap='inferno')\n",
    "# plt.xlabel(r'site $j$')\n",
    "# plt.ylabel(r'time $t$')\n",
    "# title = r'(a) QSim $\\langle N_j (t) \\rangle$: $N_s = $' + str(N_s)\n",
    "# title += r', $J = $' + str(J)\n",
    "# title += r', $h = $' + str(h)\n",
    "# title += r', $g = $' + str(g)\n",
    "# title += r', $\\Delta t=$' + str(dt)\n",
    "# plt.title(title)\n",
    "# plt.colorbar()\n",
    "# plt.savefig('qc-occ-nums.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<BarContainer object of 8 artists>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAcfElEQVR4nO3df2xV93nA4dc2tR0KOCQ0diAOXn6shKXYqV0zJ8vSrV5ZhLJm2jovyobnpkhbQaOzOrVsE+4PrWZriug6BE1akqlZBNvUpNvakjIvZIrqisQMLUnXbOlKoEltQGtt4kp2Zd/9UdWRF0i4cJ03Ns8jHQkff8+97wmx/OH4XN+yQqFQCACAJOXZAwAAFzYxAgCkEiMAQCoxAgCkEiMAQCoxAgCkEiMAQCoxAgCkmpc9wNmYnJyMF198MRYuXBhlZWXZ4wAAZ6FQKMSpU6di6dKlUV5+5usfsyJGXnzxxaivr88eAwA4B8eOHYsrrrjijJ+fFTGycOHCiPjJySxatCh5GgDgbIyMjER9ff3U9/EzmRUx8tMfzSxatEiMAMAs81q3WLiBFQBIJUYAgFRiBABIJUYAgFRiBABIJUYAgFRiBABIJUYAgFRiBABIJUYAgFRiBABIJUYAgFRiBABIJUYAgFTzsgcAuJA1fOQr2SOcsyNb12aPwBzhyggAkMqVkQvEbP3Xl395Acx9rowAAKnECACQSowAAKnECACQSowAAKnECACQSowAAKnECACQSowAAKnECACQSowAAKnECACQSowAAKnECACQ6pxiZMeOHdHQ0BDV1dWxevXqOHjw4BnX3n///VFWVjZtq66uPueBAYC5pegY2bt3b3R3d0dPT08cOnQoGhsbY82aNXH8+PEzHrNo0aL4/ve/P7U9//zz5zU0ADB3FB0j27Zti/Xr10dXV1esXLkydu3aFfPnz4/du3ef8ZiysrKoq6ub2mpra89raABg7igqRsbHx2NgYCDa29tffoDy8mhvb4/+/v4zHvfSSy/F8uXLo76+Pt7znvfEM88886rPMzY2FiMjI9M2AGBuKipGTp48GRMTE6+4slFbWxuDg4OnPeatb31r7N69O7785S/HAw88EJOTk3HjjTfG9773vTM+T29vb9TU1Ext9fX1xYwJAMwiM/5qmra2tli3bl00NTXFLbfcEl/60pfiLW95S3zuc5874zGbN2+O4eHhqe3YsWMzPSYAkGReMYuXLFkSFRUVMTQ0NG3/0NBQ1NXVndVjvOlNb4obbrghnnvuuTOuqaqqiqqqqmJGAwBmqaKujFRWVkZzc3P09fVN7ZucnIy+vr5oa2s7q8eYmJiIp556Ki6//PLiJgUA5qSiroxERHR3d0dnZ2e0tLREa2trbN++PUZHR6OrqysiItatWxfLli2L3t7eiIj4+Mc/Hj//8z8f11xzTfzwhz+MT33qU/H888/H+9///tKeCQAwKxUdIx0dHXHixInYsmVLDA4ORlNTU+zbt2/qptajR49GefnLF1x+8IMfxPr162NwcDAWL14czc3N8Y1vfCNWrlxZurMAAGatskKhUMge4rWMjIxETU1NDA8Px6JFi7LHmZUaPvKV7BHOyZGta7NHgBk1W782I3x98trO9vu396YBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKdU4zs2LEjGhoaorq6OlavXh0HDx48q+P27NkTZWVlcfvtt5/L0wIAc1DRMbJ3797o7u6Onp6eOHToUDQ2NsaaNWvi+PHjr3rckSNH4kMf+lDcfPPN5zwsADD3FB0j27Zti/Xr10dXV1esXLkydu3aFfPnz4/du3ef8ZiJiYm4884742Mf+1hcddVV5zUwADC3FBUj4+PjMTAwEO3t7S8/QHl5tLe3R39//xmP+/jHPx6XXXZZ3HXXXWf1PGNjYzEyMjJtAwDmpqJi5OTJkzExMRG1tbXT9tfW1sbg4OBpj3n88cfjC1/4Qtx7771n/Ty9vb1RU1MztdXX1xczJgAwi8zoq2lOnToVv/u7vxv33ntvLFmy5KyP27x5cwwPD09tx44dm8EpAYBM84pZvGTJkqioqIihoaFp+4eGhqKuru4V67/zne/EkSNH4rbbbpvaNzk5+ZMnnjcvnn322bj66qtfcVxVVVVUVVUVMxoAMEsVdWWksrIympubo6+vb2rf5ORk9PX1RVtb2yvWr1ixIp566qk4fPjw1PZrv/Zr8Uu/9Etx+PBhP34BAIq7MhIR0d3dHZ2dndHS0hKtra2xffv2GB0dja6uroiIWLduXSxbtix6e3ujuro6rr/++mnHX3zxxRERr9gPAFyYio6Rjo6OOHHiRGzZsiUGBwejqakp9u3bN3VT69GjR6O83C92BQDOTtExEhGxcePG2Lhx42k/d+DAgVc99v777z+XpwQA5iiXMACAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVOcUIzt27IiGhoaorq6O1atXx8GDB8+49ktf+lK0tLTExRdfHG9+85ujqakpvvjFL57zwADA3FJ0jOzduze6u7ujp6cnDh06FI2NjbFmzZo4fvz4addfcskl8ad/+qfR398f//Ef/xFdXV3R1dUVjzzyyHkPDwDMfkXHyLZt22L9+vXR1dUVK1eujF27dsX8+fNj9+7dp13/zne+M3791389rrvuurj66qtj06ZNsWrVqnj88cfPe3gAYPYrKkbGx8djYGAg2tvbX36A8vJob2+P/v7+1zy+UChEX19fPPvss/GLv/iLZ1w3NjYWIyMj0zYAYG4qKkZOnjwZExMTUVtbO21/bW1tDA4OnvG44eHhWLBgQVRWVsbatWvjs5/9bPzKr/zKGdf39vZGTU3N1FZfX1/MmADALPK6vJpm4cKFcfjw4XjiiSfiz//8z6O7uzsOHDhwxvWbN2+O4eHhqe3YsWOvx5gAQIJ5xSxesmRJVFRUxNDQ0LT9Q0NDUVdXd8bjysvL45prromIiKampvjP//zP6O3tjXe+852nXV9VVRVVVVXFjAYAzFJFXRmprKyM5ubm6Ovrm9o3OTkZfX190dbWdtaPMzk5GWNjY8U8NQAwRxV1ZSQioru7Ozo7O6OlpSVaW1tj+/btMTo6Gl1dXRERsW7duli2bFn09vZGxE/u/2hpaYmrr746xsbG4qtf/Wp88YtfjJ07d5b2TACAWanoGOno6IgTJ07Eli1bYnBwMJqammLfvn1TN7UePXo0ystfvuAyOjoaH/jAB+J73/teXHTRRbFixYp44IEHoqOjo3RnAQDMWmWFQqGQPcRrGRkZiZqamhgeHo5FixZljzMrNXzkK9kjnJMjW9dmjwAzarZ+bUb4+uS1ne33b+9NAwCkEiMAQCoxAgCkEiMAQCoxAgCkEiMAQCoxAgCkEiMAQCoxAgCkEiMAQCoxAgCkEiMAQCoxAgCkEiMAQCoxAgCkmpc9AJRSw0e+kj3COTmydW32CABpXBkBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAg1QX/0t7Z+lLQCC8HBWBucGUEAEglRgCAVGIEAEglRgCAVBf8DazAG9dsvcHczeVQHFdGAIBUrowAMONc5eLVuDICAKQSIwBAKjECAKQSIwBAKjECAKQSIwBAKjECAKQSIwBAKjECAKQSIwBAKjECAKQSIwBAKjECAKQSIwBAKjECAKQSIwBAKjECAKQSIwBAKjECAKQSIwBAKjECAKQSIwBAKjECAKQSIwBAKjECAKQ6pxjZsWNHNDQ0RHV1daxevToOHjx4xrX33ntv3HzzzbF48eJYvHhxtLe3v+p6AODCUnSM7N27N7q7u6OnpycOHToUjY2NsWbNmjh+/Php1x84cCDuuOOOePTRR6O/vz/q6+vj3e9+d7zwwgvnPTwAMPsVHSPbtm2L9evXR1dXV6xcuTJ27doV8+fPj927d592/d/+7d/GBz7wgWhqaooVK1bE5z//+ZicnIy+vr7zHh4AmP2KipHx8fEYGBiI9vb2lx+gvDza29ujv7//rB7jRz/6Ufz4xz+OSy65pLhJAYA5aV4xi0+ePBkTExNRW1s7bX9tbW18+9vfPqvH+PCHPxxLly6dFjT/39jYWIyNjU19PDIyUsyYAMAs8rq+mmbr1q2xZ8+eeOihh6K6uvqM63p7e6OmpmZqq6+vfx2nBABeT0XFyJIlS6KioiKGhoam7R8aGoq6urpXPfbuu++OrVu3xte//vVYtWrVq67dvHlzDA8PT23Hjh0rZkwAYBYpKkYqKyujubl52s2nP70Zta2t7YzH/eVf/mV84hOfiH379kVLS8trPk9VVVUsWrRo2gYAzE1F3TMSEdHd3R2dnZ3R0tISra2tsX379hgdHY2urq6IiFi3bl0sW7Ysent7IyLiL/7iL2LLli3x4IMPRkNDQwwODkZExIIFC2LBggUlPBUAYDYqOkY6OjrixIkTsWXLlhgcHIympqbYt2/f1E2tR48ejfLyly+47Ny5M8bHx+M3f/M3pz1OT09PfPSjHz2/6QGAWa/oGImI2LhxY2zcuPG0nztw4MC0j48cOXIuTwEAXCC8Nw0AkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkEqMAACpxAgAkGpe9gBA8Ro+8pXsEc7Jka1rs0cA3oBcGQEAUokRACDVOcXIjh07oqGhIaqrq2P16tVx8ODBM6595pln4jd+4zeioaEhysrKYvv27ec6KwAwBxUdI3v37o3u7u7o6emJQ4cORWNjY6xZsyaOHz9+2vU/+tGP4qqrroqtW7dGXV3deQ8MAMwtRcfItm3bYv369dHV1RUrV66MXbt2xfz582P37t2nXf+Od7wjPvWpT8Vv//ZvR1VV1XkPDADMLUXFyPj4eAwMDER7e/vLD1BeHu3t7dHf31+yocbGxmJkZGTaBgDMTUXFyMmTJ2NiYiJqa2un7a+trY3BwcGSDdXb2xs1NTVTW319fckeGwB4Y3lDvppm8+bNMTw8PLUdO3YseyQAYIYU9UvPlixZEhUVFTE0NDRt/9DQUElvTq2qqnJ/CQBcIIq6MlJZWRnNzc3R19c3tW9ycjL6+vqira2t5MMBAHNf0b8Ovru7Ozo7O6OlpSVaW1tj+/btMTo6Gl1dXRERsW7duli2bFn09vZGxE9uev3Wt7419ecXXnghDh8+HAsWLIhrrrmmhKcCAMxGRcdIR0dHnDhxIrZs2RKDg4PR1NQU+/btm7qp9ejRo1Fe/vIFlxdffDFuuOGGqY/vvvvuuPvuu+OWW26JAwcOnP8ZAACz2jm9Ud7GjRtj48aNp/3c/w+MhoaGKBQK5/I0AMAF4A35ahoA4MIhRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVGIEAEglRgCAVPOyBwCAuaLhI1/JHuGcHNm6NvX5XRkBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAglRgBAFKJEQAg1TnFyI4dO6KhoSGqq6tj9erVcfDgwVdd//d///exYsWKqK6ujre97W3x1a9+9ZyGBQDmnqJjZO/evdHd3R09PT1x6NChaGxsjDVr1sTx48dPu/4b3/hG3HHHHXHXXXfFv//7v8ftt98et99+ezz99NPnPTwAMPsVHSPbtm2L9evXR1dXV6xcuTJ27doV8+fPj927d592/Wc+85n41V/91fjjP/7juO666+ITn/hEvP3tb4+//uu/Pu/hAYDZb14xi8fHx2NgYCA2b948ta+8vDza29ujv7//tMf09/dHd3f3tH1r1qyJhx9++IzPMzY2FmNjY1MfDw8PR0TEyMhIMeOelcmxH5X8MV8vxfz3mK3nWezfufN8Y3OerzRbzzHiwjhP/8+W5nELhcKrrisqRk6ePBkTExNRW1s7bX9tbW18+9vfPu0xg4ODp10/ODh4xufp7e2Nj33sY6/YX19fX8y4c17N9uwJZt6FcI4RznOucZ5zx4VwjhEzf56nTp2KmpqaM36+qBh5vWzevHna1ZTJycn43//937j00kujrKwscbKzNzIyEvX19XHs2LFYtGhR9jgzxnnOLc5z7rgQzjHCeb7RFQqFOHXqVCxduvRV1xUVI0uWLImKiooYGhqatn9oaCjq6upOe0xdXV1R6yMiqqqqoqqqatq+iy++uJhR3zAWLVo0q/7HOVfOc25xnnPHhXCOEc7zjezVroj8VFE3sFZWVkZzc3P09fVN7ZucnIy+vr5oa2s77TFtbW3T1kdE7N+//4zrAYALS9E/punu7o7Ozs5oaWmJ1tbW2L59e4yOjkZXV1dERKxbty6WLVsWvb29ERGxadOmuOWWW+LTn/50rF27Nvbs2RNPPvlk3HPPPaU9EwBgVio6Rjo6OuLEiROxZcuWGBwcjKampti3b9/UTapHjx6N8vKXL7jceOON8eCDD8af/dmfxZ/8yZ/EtddeGw8//HBcf/31pTuLN6Cqqqro6el5xY+b5hrnObc4z7njQjjHCOc5V5QVXuv1NgAAM8h70wAAqcQIAJBKjAAAqcQIAJBKjMyQHTt2RENDQ1RXV8fq1avj4MGD2SOV1L/927/FbbfdFkuXLo2ysrJXfa+h2ay3tzfe8Y53xMKFC+Oyyy6L22+/PZ599tnssUpu586dsWrVqqlfqNTW1hZf+9rXsseaUVu3bo2ysrL44Ac/mD1KSX30ox+NsrKyaduKFSuyx5oRL7zwQvzO7/xOXHrppXHRRRfF2972tnjyySezxyqphoaGV/x9lpWVxYYNG7JHKykxMgP27t0b3d3d0dPTE4cOHYrGxsZYs2ZNHD9+PHu0khkdHY3GxsbYsWNH9igz6rHHHosNGzbEN7/5zdi/f3/8+Mc/jne/+90xOjqaPVpJXXHFFbF169YYGBiIJ598Mn75l3853vOe98QzzzyTPdqMeOKJJ+Jzn/tcrFq1KnuUGfFzP/dz8f3vf39qe/zxx7NHKrkf/OAHcdNNN8Wb3vSm+NrXvhbf+ta34tOf/nQsXrw4e7SSeuKJJ6b9Xe7fvz8iIt773vcmT1ZiBUqutbW1sGHDhqmPJyYmCkuXLi309vYmTjVzIqLw0EMPZY/xujh+/HghIgqPPfZY9igzbvHixYXPf/7z2WOU3KlTpwrXXnttYf/+/YVbbrmlsGnTpuyRSqqnp6fQ2NiYPcaM+/CHP1z4hV/4hewxXnebNm0qXH311YXJycnsUUrKlZESGx8fj4GBgWhvb5/aV15eHu3t7dHf3584GaUwPDwcERGXXHJJ8iQzZ2JiIvbs2ROjo6Nz8m0bNmzYEGvXrp32NTrX/Pd//3csXbo0rrrqqrjzzjvj6NGj2SOV3D/+4z9GS0tLvPe9743LLrssbrjhhrj33nuzx5pR4+Pj8cADD8T73ve+WfOmsWdLjJTYyZMnY2JiYuo30v5UbW1tDA4OJk1FKUxOTsYHP/jBuOmmm+bkbxB+6qmnYsGCBVFVVRW///u/Hw899FCsXLkye6yS2rNnTxw6dGjq7SrmotWrV8f9998f+/bti507d8Z3v/vduPnmm+PUqVPZo5XU//zP/8TOnTvj2muvjUceeST+4A/+IP7wD/8w/uZv/iZ7tBnz8MMPxw9/+MP4vd/7vexRSq7oXwcPF6oNGzbE008/PSd//h4R8da3vjUOHz4cw8PD8Q//8A/R2dkZjz322JwJkmPHjsWmTZti//79UV1dnT3OjLn11lun/rxq1apYvXp1LF++PP7u7/4u7rrrrsTJSmtycjJaWlrik5/8ZERE3HDDDfH000/Hrl27orOzM3m6mfGFL3whbr311li6dGn2KCXnykiJLVmyJCoqKmJoaGja/qGhoairq0uaivO1cePG+Od//ud49NFH44orrsgeZ0ZUVlbGNddcE83NzdHb2xuNjY3xmc98JnuskhkYGIjjx4/H29/+9pg3b17MmzcvHnvssfirv/qrmDdvXkxMTGSPOCMuvvji+Nmf/dl47rnnskcpqcsvv/wVoXzdddfNyR9JRUQ8//zz8S//8i/x/ve/P3uUGSFGSqyysjKam5ujr69vat/k5GT09fXNyZ+/z3WFQiE2btwYDz30UPzrv/5r/MzP/Ez2SK+bycnJGBsbyx6jZN71rnfFU089FYcPH57aWlpa4s4774zDhw9HRUVF9ogz4qWXXorvfOc7cfnll2ePUlI33XTTK15m/1//9V+xfPnypIlm1n333ReXXXZZrF27NnuUGeHHNDOgu7s7Ojs7o6WlJVpbW2P79u0xOjoaXV1d2aOVzEsvvTTtX1rf/e534/Dhw3HJJZfElVdemThZaW3YsCEefPDB+PKXvxwLFy6cuu+npqYmLrroouTpSmfz5s1x6623xpVXXhmnTp2KBx98MA4cOBCPPPJI9mgls3Dhwlfc6/PmN785Lr300jl1D9CHPvShuO2222L58uXx4osvRk9PT1RUVMQdd9yRPVpJ/dEf/VHceOON8clPfjJ+67d+Kw4ePBj33HNP3HPPPdmjldzk5GTcd9990dnZGfPmzdFv29kv55mrPvvZzxauvPLKQmVlZaG1tbXwzW9+M3ukknr00UcLEfGKrbOzM3u0kjrdOUZE4b777sseraTe9773FZYvX16orKwsvOUtbym8613vKnz961/PHmvGzcWX9nZ0dBQuv/zyQmVlZWHZsmWFjo6OwnPPPZc91oz4p3/6p8L1119fqKqqKqxYsaJwzz33ZI80Ix555JFCRBSeffbZ7FFmTFmhUCjkZBAAgHtGAIBkYgQASCVGAIBUYgQASCVGAIBUYgQASCVGAIBUYgQASCVGAIBUYgQASCVGAIBUYgQASPV/GeFS376si/sAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(range(N_s), qc_occ_nums[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Larger systems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [],
   "source": [
    "exact_occs = np.loadtxt(\"ising-approx-occs-N_s12-J0.4-h1-g0.3-width1.5-dt0.6.txt\", dtype=float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(15, 12)"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exact_occs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'Y' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[113], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m plt\u001b[38;5;241m.\u001b[39mpcolormesh(X, \u001b[43mY\u001b[49m, exact_occs, cmap\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124minferno\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m      2\u001b[0m plt\u001b[38;5;241m.\u001b[39mxlabel(\u001b[38;5;124mr\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124msite $j$\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m      3\u001b[0m plt\u001b[38;5;241m.\u001b[39mylabel(\u001b[38;5;124mr\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mtime $t$\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "\u001b[0;31mNameError\u001b[0m: name 'Y' is not defined"
     ]
    }
   ],
   "source": [
    "plt.pcolormesh(X, Y, exact_occs, cmap='inferno')\n",
    "plt.xlabel(r'site $j$')\n",
    "plt.ylabel(r'time $t$')\n",
    "title = r'(a) QSim $\\langle N_j (t) \\rangle$: $N_s = $' + str(N_s)\n",
    "title += r', $J = $' + str(J)\n",
    "title += r', $h = $' + str(h)\n",
    "title += r', $g = $' + str(g)\n",
    "title += r', $\\Delta t=$' + str(dt)\n",
    "plt.title(title)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'Y' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[114], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m plt\u001b[38;5;241m.\u001b[39mpcolormesh(X, \u001b[43mY\u001b[49m, np\u001b[38;5;241m.\u001b[39mabs(exact_occs \u001b[38;5;241m-\u001b[39m qc_occ_nums), cmap\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124minferno\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m      2\u001b[0m plt\u001b[38;5;241m.\u001b[39mxlabel(\u001b[38;5;124mr\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124msite $j$\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m      3\u001b[0m plt\u001b[38;5;241m.\u001b[39mylabel(\u001b[38;5;124mr\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mtime $t$\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "\u001b[0;31mNameError\u001b[0m: name 'Y' is not defined"
     ]
    }
   ],
   "source": [
    "plt.pcolormesh(X, Y, np.abs(exact_occs - qc_occ_nums), cmap='inferno')\n",
    "plt.xlabel(r'site $j$')\n",
    "plt.ylabel(r'time $t$')\n",
    "title = r'(a) QSim $\\langle N_j (t) \\rangle$: $N_s = $' + str(N_s)\n",
    "title += r', $J = $' + str(J)\n",
    "title += r', $h = $' + str(h)\n",
    "title += r', $g = $' + str(g)\n",
    "title += r', $\\Delta t=$' + str(dt)\n",
    "plt.title(title)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 341,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.13212111247685732"
      ]
     },
     "execution_count": 341,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(np.abs(exact_occs - qc_occ_nums))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for i in range(N_t):\n",
    "#     plt.bar(range(N_s), qc_occ_nums[i], label='Exact Trotter', color='gray')\n",
    "#     plt.xlabel('s (site)')\n",
    "#     plt.ylabel(r'$\\langle n_s \\rangle$')\n",
    "#     plt.legend()\n",
    "#     plt.ylim(0,1.3)\n",
    "#     title = r'$N_s =$' + str(N_s)\n",
    "#     title += r', $\\lambda$ = ' + str(lam)\n",
    "#     title += r', $t = $' + str(round(i * dt, 4))\n",
    "#     plt.title(title)\n",
    "#     folder = 'qc-scattering-plots/'\n",
    "#     filename = folder\n",
    "#     filename += 'scatt-timestep' + str(i) + '.png'\n",
    "#     plt.savefig(filename)\n",
    "#     plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exact Diagonalization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "CDAG = []\n",
    "for s in range(N_s):\n",
    "    cdag = Sm[s]\n",
    "    for j in range(s):\n",
    "        cdag = -Z[j] @ cdag\n",
    "    CDAG.append(cdag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([-1, 0, 1, 2])\n",
    "def GWP_coeff(x, x_mean, x_dif, nk_mean, sig_k):\n",
    "    k = nk_mean * np.pi / N_s\n",
    "    coeff = 1 / np.sqrt(N_s) * np.exp(-1j * k * (x + x_mean))\n",
    "    coeff *= np.exp(-x_dif[x]**2 / sig_k**2)\n",
    "    return coeff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "xr_mean = 1\n",
    "xl_mean = 5\n",
    "nk_mean = 3\n",
    "sig_k = 3/2\n",
    "rm_op = np.zeros((2**N_s, 2**N_s), dtype=complex)\n",
    "lm_op = np.zeros((2**N_s, 2**N_s), dtype=complex)\n",
    "for j in range(4):\n",
    "    rm_op += GWP_coeff(j, 1, x, nk_mean, sig_k) * CDAG[j]\n",
    "    lm_op += GWP_coeff(j, 5, x, -nk_mean, sig_k) * CDAG[j+4]\n",
    "\n",
    "exact_init_state = np.zeros((2**N_s), dtype=complex)\n",
    "exact_init_state[0] = 1\n",
    "exact_init_state = lm_op @ rm_op @ ground_state\n",
    "exact_init_state /= np.linalg.norm(exact_init_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "state = np.array(exact_init_state)\n",
    "occs = []\n",
    "for t in range(N_t):\n",
    "    occ = []\n",
    "    for j in range(N_s):\n",
    "        occ_j = np.inner(state.conjugate(), N[j] @ state)\n",
    "        occ.append(occ_j)\n",
    "    occs.append(np.real(occ))\n",
    "    state = TEO @ state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Initial State: $N_s = $8, $n_k = \\\\pm$3, $\\\\sigma_k=$1.5')"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAHNCAYAAAAaKaG7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABAzklEQVR4nO3de1yUZf7/8fcwnMSziOApyNQ8o4IammkbSpaalqZ+LYnKdjcpi2r7Wa1olriVZQdX01LbykOlttV6jM3agjwReD5rmgqIJxQVcOb+/dE62ygY6gwz3Lyej8f9eMQ1133P5x5j5s11X9c9FsMwDAEAAJiYj6cLAAAAcDcCDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDyq0OXPmyGKxaN++fS7t68p9AQCeR+CBx10IE+vWrXPJ8dLS0jRu3DidOHHCJce7Uhs3btSgQYMUHh6uwMBANWzYUL169dLbb7/t0ho9cZ7Lly+XxWKRxWLRli1bLnm8X79+atSoUbnV81s7d+7U0KFD1ahRIwUFBalFixZ68cUXdebMGY/UU5Fs3rxZgwcPVpMmTRQUFKS6devqlltu0Zdffunp0lzu9OnTSk5O1u233646derIYrFozpw5Zd5/1apVjt+Bi7cff/zRfYXjmhF4UKHdf//9Onv2rMLDwx1taWlpGj9+/CVBoKS+rpaWlqbo6GhlZWVp5MiReuedd/Twww/Lx8dHb7755u/WeKXPda3HuFJZWVmSJB8fH3311VclPt6uXbtyq+eCAwcOqHPnzvrxxx+VmJioKVOmKCYmRsnJyRo2bFi511PR/Pzzzzp16pTi4+P15ptv6q9//askqX///poxY4aHq3OtvLw8vfjii9q6dasiIyOv+jiPP/64PvzwQ6etadOmLqwUrubr6QKAa2G1WmW1Wl3e92q9/PLLqlmzptauXatatWo5PZabm+vW5y4PGzZsUI0aNdStWzd9+eWX+stf/uJ47Pjx4zpw4ICGDx9e7nV9+OGHOnHihL7//nu1bt1akvTII4/IbrfrH//4h44fP67atWuXe12e0rNnT0VERJR55OKOO+7QHXfc4dSWmJioqKgovf7663rkkUfcUKVn1K9fX4cPH1ZYWJjWrVunTp06XdVxunfvrkGDBrm4OrgTIzzwOuPGjZPFYtGuXbv0wAMPqFatWqpZs6YSEhIuuTxx8dyacePG6ZlnnpEkXX/99Y6h5n379pU4D+fnn3/Wo48+qhtvvFFVqlRRcHCwBg8efNVzdXbv3q3WrVtfEnYkqV69er9bY1lr+r1jSNLBgwf14IMPKjQ0VAEBAWrdurVmzZp1SV3btm3T/v37y3R+WVlZatu2rfr27av09HQdPXrU6TFJHhnhyc/PlySFhoY6tdevX18+Pj7y9/e/quPeeuutuuWWW5SRkaE+ffqoevXqatiwodNonVlZrVY1btz4mkcQZ8yYoY4dOyooKOiSS0BNmjRxTbFXICAgQGFhYS451qlTp3T+/HmXHAvuxwgPvNa9996r66+/XikpKcrIyNB7772nevXq6W9/+1up+9x9993asWOH5s2bpzfeeEN169aVJIWEhJTYf+3atUpLS3PM/di3b5+mTZumnj17asuWLQoKCrqimsPDw5Wenq5NmzapTZs2V1VjWWr6vWPk5OTopptuksViUWJiokJCQrR06VI99NBDys/P1xNPPOGop2XLlurRo4dWrVp12XMrKirS9u3bNXLkSPXt21ejRo3SkiVLdP/990v6dfRH0hVfJiguLtbJkyfL1LdOnTry8bn077SePXvqb3/7mx566CGNHz9ewcHBSktL07Rp0/T444+ratWqV1TTBRs3blSDBg3Ur18/JSQkaMCAAZo5c6aefPJJ/eEPf1Dbtm2v6HiuONeyHLO4uFiFhYXKy8u7omMWFBTo7NmzOnnypL744gstXbpUQ4YMKVO9JXnyySc1ZcoU9e7dWwkJCfrll1/0xhtvqLi4WH379lVUVNQVHc8dr9/VSkhI0OnTp2W1WtW9e3e9+uqrio6OdtvzwQUMwMNmz55tSDLWrl1rGIZhJCcnG5KMBx980KnfwIEDjeDg4BL33bt3r6Pt1VdfvaSttL5nzpy5pJ709HRDkvGPf/zjsvuWZMWKFYbVajWsVqsRExNj/OUvfzGWL19uFBUVOfUrrcYrqelyx3jooYeM+vXrG3l5eU7tQ4cONWrWrOn0HJKMHj16XPa8DMMwfvrpJ0OSMX36dMMwDKNt27bG4MGDHY8/+OCDRkBAgHH+/PnfPdZvffPNN4akMm2Xe/0nTJhgVKlSxan/888/f0W1/NahQ4cMSUZISIhx4MABR/uWLVsMScYHH3xwxcd01bm665h//OMfHX19fHyMQYMGGceOHbvi8zQMw/juu+8MScaf//xnp/bx48cbkow1a9Zc8TFd/fqtXbvWkGTMnj27zDX88MMPxj333GO8//77xj//+U8jJSXFCA4ONgIDA42MjIwrPieUH0Z44LX+9Kc/Of3cvXt3LV68WPn5+apRo4ZLnqNKlSqO/y4uLlZ+fr6aNm2qWrVqKSMjwzF6UVa9evVSenq6UlJStHz5cqWnp+uVV15RSEiI3nvvPfXv39/tNRmGoYULF+ree++VYRhOf+XHxcVp/vz5ysjIULdu3Rz9y+LCCM6FS1Z9+/bV1KlTVVxcLD8/P2VlZal169ZXPE8qMjJSK1euLFPfy12KiIiI0C233KJ77rlHwcHB+te//qWJEycqLCxMiYmJV1ST9OvojiQlJyc7rTzz8/OTJKfLZIZhqHr16tqzZ4/j0mVJXHWuv3fMp556SmFhYY7LnmU95hNPPKFBgwbp0KFD+uSTT2Sz2VRUVFSmOi72xhtvqE6dOnr11Ved2nv06CFJ2rFjh2P+jCdfvyvVtWtXde3a1fFz//79NWjQILVr105jxozRsmXL3PK8uHYEHnit6667zunnC5NOjx8/7rLAc/bsWaWkpGj27Nk6ePCg04d/WYfOL9apUyctWrRIRUVFysrK0uLFi/XGG29o0KBByszMVKtWrdxa05EjR3TixAnNmDGj1BU2VzOBOisrSxaLxXEZp2/fvkpJSdF3332nnj17avPmzRo6dOgVH7d27dqKjY294v1+a/78+XrkkUe0Y8cORzi5++67Zbfb9eyzz2rYsGEKDg6+omNeCDwDBgxwat+2bZsk6cYbb3S07d27V0FBQZf9sJZcc65lOWbt2rVVv379K36uFi1aqEWLFpKkESNGqHfv3urXr59Wr14ti8VS5uOcP39eK1eu1F133XXJ5cQLAeq3v8OefP1coWnTprrrrru0aNEi2Ww2ty+OwNUh8MBrlfamUdYRibJ47LHHNHv2bD3xxBOKiYlRzZo1ZbFYNHToUNnt9ms6tr+/vzp16qROnTqpefPmSkhI0Keffqrk5GS31nShz3333af4+PgS+1zNxOINGzaoSZMmqlatmiTppptuUt26dfXll1+qQYMGOnfu3FUt8y0qKtKxY8fK1DckJKTE/y/+/ve/q0OHDpfcA6h///6aM2eOfvrppyv+oNywYYPCwsLUsGFDp/asrCz5+vo6BddNmzY5VoddjivOtTwNGjRIf/zjH7Vjxw6ngPd79u3bp9OnT5c4j239+vWSfp07doEZXr/GjRurqKhIBQUFLvuDDK5F4IHpXMlfop999pni4+M1efJkR9u5c+dcfm+bC5MZDx8+/Ls1lrWm0o4REhKi6tWry2azufSv4Q0bNjgug0m/3ounT58++vLLL3XTTTdJcg5S27ZtU2JiojIyMmQYhu6//3699dZblxw3LS1Nt956a5lq2Lt3ryIiIi5pz8nJKXHZeXFxsSRd1UqajRs3lhjgNmzYoObNmysgIMDR9tsP7NzcXN1zzz3q2bOnXnzxRad/J1eca3k6e/aspCsf7Tx16pQkXbI6zjAMffrpp2rdurXTPWvM8Prt2bNHgYGBjj8I4H0IPDCdC0PoZQktVqv1khGjt99+Wzab7aqe+5tvvlHPnj0vCSNLliyR9L/LIJersaw1lXYMq9Wqe+65R3Pnzi1xtdiRI0ecVq1t27ZNQUFBl1xC/K3s7Gzl5uZeEgD69u2rDz/8UPPmzZPkvEJr+PDhevbZZzV48GCdOnVKO3fuLPHYrpiX0bx5c61YsUI7duxQ8+bNHe3z5s2Tj4/PFY9o2Ww2bd26Vb169brksaysLHXo0MGpbdOmTerevbt++uknDRkyROPGjdP//d//XbJvec1B+b0VdxfLzc295HJScXGx/vGPf6hKlSq/exn2Yhf+X/r666+VlJTkaJ8yZYoyMjL00UcfOfX3ttfvgjNnzmj//v2qW7euYyXkxb8/0q//T3zxxRfq06ePW1eG4doQeGA6F5a6Pv/88xo6dKj8/PzUr1+/Evte+MCuWbOmWrVqpfT0dH399ddXPN/jgscee0xnzpzRwIED1aJFCxUVFSktLU0LFixQRESEEhISLltj1apVy1zT5Y4xadIkffPNN+rSpYtGjhypVq1a6dixY8rIyNDXX3/tdFmgLMvSS7vHTlxcnPz8/ByXtX5b4+7du1VUVCS73a4aNWqUugTZFfMynnnmGS1dulTdu3dXYmKigoOD9dVXX2np0qV6+OGH1aBBA6f+Fovlsue8c+fOEi/RnT17Vrt27brkUuGmTZscE3Tnzp3rGPG6mDvmoOTk5JQ5BAwcOLDEJfp//OMflZ+fr1tuuUUNGzZUdna2Pv74Y23btk2TJ092GrX4vddOkoKDgzVgwAB9/vnnGj58uLp166bvv/9e8+bN08MPP3zJzSnL+/V75513dOLECR06dEiS9OWXX+qXX36R9OvvcM2aNSVJa9as0a233qrk5GSNGzdOkjRkyBBVqVJFXbt2Vb169bRlyxbNmDFDQUFBmjRp0jXXBjfyzOIw4H9KW5Z+5MiREvv9drlpacvFJ0yYYDRs2NDw8fFxPF5S3+PHjxsJCQlG3bp1jWrVqhlxcXHGtm3bjPDwcCM+Pv53n+diS5cuNR588EGjRYsWRrVq1Qx/f3+jadOmxmOPPWbk5OT8bo1XUtPljmEYhpGTk2OMGjXKaNy4seHn52eEhYUZt912mzFjxgynY6gMy9JfeeUVQ5Kxa9euSx679dZbDUnG7bff7tS+ZMkSo1u3bkZoaKjxzDPPGMXFxZd9jmu1evVqo0+fPkZYWJjh5+dnNG/e3Hj55Zcved5Tp04ZkoyhQ4eWeqxPPvnEkGRs2rTJqX3NmjWGJOOrr75ytBUXFxv+/v5GSEiI020DyosrlmrPmzfPiI2NNUJDQw1fX1+jdu3aRmxsrPHPf/7TqV9ZXrsLjh8/bjzwwANG7dq1jYCAAKNDhw7G+++/f0k/T7x+4eHhZXqNLry2ycnJjrY333zT6Ny5s1GnTh3D19fXqF+/vnHfffcZO3fuLJfacfUshuHCGaAAcJGff/5Zt9xyi2bOnKnevXt7uhwtWbJEffv2ddw1+lpt2bJFvXr10ieffKLBgwcrPT3drd/X5kmufu2kyvX6wbO42AjA5RYuXKi9e/dK+vU2AkVFRY7lzp72zTffaOjQoS77wN60aZPatm2rbt266bnnntPdd9+tc+fOueTY3sbVr51UuV4/eBYjPABc7vHHH9cnn3yi06dP64YbbtDEiRN15513erostxg7dqzOnTunV155RZIc83s++OADT5ZVYfD6obwQeAAAgOlxSQsAAJgegQcAAJgegQcAAJgeNx78L7vdrkOHDql69epX9NUEAADAcwzD0KlTp9SgQYPL3umawPNfhw4dUuPGjT1dBgAAuAoHDhy45AuEf4vA81/Vq1eX9OsLxjfdAgBQMeTn56tx48aOz/HSEHj+68JlrBo1ahB4AACoYH5vOgqTlgEAgOkReAAAgOkReAAAgOkxhwcAABex2WwqLi72dBmm4ufnJ6vVes3HIfAAAHCNDMNQdna2Tpw44elSTKlWrVoKCwu7pvvkEXgAALhGF8JOvXr1FBQUxA1sXcQwDJ05c0a5ubmSpPr161/1sQg8AABcA5vN5gg7wcHBni7HdKpUqSJJys3NVb169a768haTlgEAuAYX5uwEBQV5uBLzuvDaXsv8KAIPAAAuwGUs93HFa0vgAQAApkfgAQDAi5y32T1dgikReAAA8BI/7jmqDhNWavWeo+X2nAcOHNCDDz6oBg0ayN/fX+Hh4Ro9erSOHnWuYdeuXUpISFCjRo0UEBCg66+/XsOGDdO6devKrdZrQeABAMALFNvsGrNoo06dO68xizequBxGevbs2aPo6Gjt3LlT8+bN065duzR9+nSlpqYqJiZGx44dkyStW7dOUVFR2rFjh959911t2bJFixcvVosWLfTUU0+5vU5XYFk6AABe4IO0fdqbVyBJ2nOkQB+k7dPD3Zu49TlHjRolf39/rVixwrH8+7rrrlOHDh10ww036Pnnn9ff//53PfDAA2rWrJn+85//yMfnf2Ml7du31+jRo91ao6swwgMAgIfl5J/Tayu2O7VNXrFDOfnn3Pacx44d0/Lly/Xoo486ws4FYWFhGj58uBYsWKDMzExt3rxZTz31lFPYuaBWrVpuq9GVGOGBy4wfP97TJVyV5ORkT5cAuJU7fzfthuTjptXYlel386WvtqjYZji1FdnsevlfW/XWsA5uec6dO3fKMAy1bNmyxMdbtmyp48ePa+fOnZKkFi1auKWO8sIIDyoEu/H7fQCUr2xbNc09117ZtmqeLqVCS9udpy83HJbtojc6m93QF1mHlL7bvROYDePyb7C/93hFQeCB1+NNFfA+dsOiH4rDVSxfpRWHy25w072rtSjjYKmjZD4WaWHGL2553qZNm8pisWjr1q0lPr5161aFhISoefPmkqRt27a5pY7yQuCBV+NNFfBOW87XU74RKEk6aQRqq62ehyuquO7p2KjUUWy7IQ2KauSW5w0ODlavXr3097//XWfPnnV6LDs7Wx9//LEeeOABtW/fXq1atdLkyZNlt1+6cqyifEM8gQdejTdVwPucMfyUcb6BpAt/gFiUUdxAZww/T5ZVYcXcEKx+7erLetEwj9XHov6RDXRTE/d9Iek777yjwsJCxcXF6bvvvtOBAwe0bNky9erVS82bN9fYsWNlsVg0e/Zs7dixQ927d9eSJUu0Z88ebdiwQS+//LLuuusut9XnSgQeeC3eVAHvtKaokewXfXzY5KM1Re4ZiagMXujbSn5W58Djb/XR83eWPKHYVZo1a6a1a9eqSZMmuvfeexUeHq4+ffqoefPm+uGHH1St2q9TCTp37qx169apadOmGjlypFq2bKn+/ftr8+bNmjJliltrdBUCD7wWb6qA9zlsq6699mAZcv5wNmTRXnuwDtuqe6iyii20RqCe7n2jU9tTvZsrtEag2587IiJCc+bMUXZ2tux2u8aOHasVK1Zow4YNTv2aN2+uDz74QAcPHlRhYaH27dunuXPnqkMH96wiczUCD7wSb6qAd9plC5ZU2qod47+P42rEd41Qk5CqkqQmIVUV3zXCI3WMHz9eb731ln788ccS5+xUVNyHB17pf2+qJU1S/vVNtb71VDlXBaCp9ah22eqW8qhFzax55VqPmfhZfTRxYFuN/Mc6pQxsKz+r58YkEhISPPbc7sIID7xSU+tRlRx2JN5UAc+pbz2l632OynLRKI9Fhq73Oaow62kPVWYONzUJ1k9/7aUubpyoXFkReOCVeFMFvFdn/1/kI+dLHVbZ1dnfPfeLqWx8PTiyY2a8qvBavKkC3inIUqyOvof0v7k8hjr6HVKQpdiTZQGXReCB1+JNFfBerXxzVcPy6xdb1rScU0trrocrAi6PwAOvxpsq4J18LIa6+f0sP51XV7+f5WMxx/ctwbwIPPBqvKkC3ivMelr/F5jJnDpUCCxLh9e78KZa2pfrAfAcfi9RUTDCgwqBN1UAwLVghAcAADcYP358uT5fcnLyFe/zwAMP6IMPPrikPS4uTsuWLXNFWZc1btw4ff7558rMzHT7cxF4AACoxG6//XbNnj3bqS0gIMBD1bgPl7QAAKjEAgICFBYW5rTVrl1bq1atkr+/v/7zn/84+r7yyiuqV6+ecnJyJEnLli3TzTffrFq1aik4OFh9+/bV7t27nY7/yy+/aNiwYapTp46qVq2q6OhorV69WnPmzNH48eOVlZUli8Uii8WiOXPmuO08GeEBAACX6Nmzp5544gndf//9ysrK0p49e/TXv/5Vn376qUJDQyVJBQUFSkpKUrt27XT69GmNHTtWAwcOVGZmpnx8fHT69Gn16NFDDRs21BdffKGwsDBlZGTIbrdryJAh2rRpk5YtW6avv/5aklSzZk23nQ+BBwCASuyrr75StWrVnNqee+45Pffcc3rppZe0cuVKPfLII9q0aZPi4+PVv39/R7977rnHab9Zs2YpJCREW7ZsUZs2bTR37lwdOXJEa9euVZ06dSRJTZs2dfSvVq2afH19FRYW5sYz/BWBBwCASuzWW2/VtGnTnNouhBN/f399/PHHateuncLDw/XGG2849du5c6fGjh2r1atXKy8vT3b7r18HtH//frVp00aZmZnq0KGD43ie5LVzeKZOnaqIiAgFBgaqS5cuWrNmTal958yZ47j+d2ELDAwsx2oBAKiYqlatqqZNmzptvw0oaWlpkqRjx47p2LFjTvv269dPx44d08yZM7V69WqtXr1aklRUVCRJqlKlSjmdxe/zysCzYMECJSUlKTk5WRkZGYqMjFRcXJxyc0v/WoEaNWro8OHDju3nn38ux4oBADCf3bt368knn9TMmTPVpUsXxcfHO0Zxjh49qu3bt+uFF17QbbfdppYtW+r48eNO+7dr106ZmZmXBKUL/P39ZbPZ3H4ekpcGntdff10jR45UQkKCWrVqpenTpysoKEizZs0qdR+LxeI0w/zChCoAAFC6wsJCZWdnO215eXmy2Wy67777FBcXp4SEBM2ePVsbNmzQ5MmTJUm1a9dWcHCwZsyYoV27dunf//63kpKSnI49bNgwhYWFacCAAfrhhx+0Z88eLVy4UOnp6ZKkiIgI7d27V5mZmcrLy1NhYaHbztPrAk9RUZHWr1+v2NhYR5uPj49iY2MdL1BJTp8+rfDwcDVu3Fh33XWXNm/efNnnKSwsVH5+vtMGAEBls2zZMtWvX99pu/nmm/Xyyy/r559/1rvvvitJql+/vmbMmKEXXnhBWVlZ8vHx0fz587V+/Xq1adNGTz75pF599VWnY/v7+2vFihWqV6+e7rjjDrVt21aTJk2S1WqV9Ouk59tvv1233nqrQkJCNG/ePLedp9dNWr6QKi8eoQkNDdW2bdtK3OfGG2/UrFmz1K5dO508eVKvvfaaunbtqs2bN6tRo0Yl7pOSklLud8EEAFQeV3Pn4/I2Z86cy977ZuzYsU4/33333U6jMLGxsdqyZYtTH8Nw/pLn8PBwffbZZyUePyAgoNTHXM3rRniuRkxMjEaMGKH27durR48eWrRokUJCQhyptCRjxozRyZMnHduBAwfKsWIAAFCevG6Ep27durJarY67OF6Qk5NT5nX6fn5+6tChg3bt2lVqn4CAAFPeOhsAAFzK60Z4/P39FRUVpdTUVEeb3W5XamqqYmJiynQMm82mjRs3qn79+u4qEwAAVCBeN8IjSUlJSYqPj1d0dLQ6d+6sKVOmqKCgQAkJCZKkESNGqGHDhkpJSZEkvfjii7rpppvUtGlTnThxQq+++qp+/vlnPfzww548DQAA4CW8MvAMGTJER44c0dixY5Wdna327dtr2bJljonM+/fvl4/P/wanjh8/rpEjRyo7O1u1a9dWVFSU0tLS1KpVK0+dAgCgkrl4si5cxxWvrVcGHklKTExUYmJiiY+tWrXK6ec33njjkttdAwBQHvz8/CRJZ86c8ao7C5vJmTNnJP3vtb4aXht4AACoCKxWq2rVquX4NoCgoCBZLBYPV2UOhmHozJkzys3NVa1atRz377kaBB4AAK7RhVXEl/sKJFy9WrVqXfM3qhN4AAC4RhaLRfXr11e9evVUXFzs6XJMxc/P75pGdi4g8AAA4CJWq9UlH85wPa+7Dw8AAICrEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpEXgAAIDpeW3gmTp1qiIiIhQYGKguXbpozZo1Zdpv/vz5slgsGjBggHsLBAAAFYZXBp4FCxYoKSlJycnJysjIUGRkpOLi4pSbm3vZ/fbt26enn35a3bt3L6dKAQBAReCVgef111/XyJEjlZCQoFatWmn69OkKCgrSrFmzSt3HZrNp+PDhGj9+vJo0aVKO1QIAAG/ndYGnqKhI69evV2xsrKPNx8dHsbGxSk9PL3W/F198UfXq1dNDDz1UpucpLCxUfn6+0wYAAMzJ6wJPXl6ebDabQkNDndpDQ0OVnZ1d4j7ff/+93n//fc2cObPMz5OSkqKaNWs6tsaNG19T3QAAwHt5XeC5UqdOndL999+vmTNnqm7dumXeb8yYMTp58qRjO3DggBurBAAAnuTr6QIuVrduXVmtVuXk5Di15+TkKCws7JL+u3fv1r59+9SvXz9Hm91ulyT5+vpq+/btuuGGGy7ZLyAgQAEBAS6uHgAAeCOvG+Hx9/dXVFSUUlNTHW12u12pqamKiYm5pH+LFi20ceNGZWZmOrb+/fvr1ltvVWZmJpeqAACA943wSFJSUpLi4+MVHR2tzp07a8qUKSooKFBCQoIkacSIEWrYsKFSUlIUGBioNm3aOO1fq1YtSbqkHQAAVE5eGXiGDBmiI0eOaOzYscrOzlb79u21bNkyx0Tm/fv3y8fH6wanAACAl/LKwCNJiYmJSkxMLPGxVatWXXbfOXPmuL4gAABQYTFMAgAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATI/AAwAATM9rA8/UqVMVERGhwMBAdenSRWvWrCm176JFixQdHa1atWqpatWqat++vT788MNyrBYAAHgzrww8CxYsUFJSkpKTk5WRkaHIyEjFxcUpNze3xP516tTR888/r/T0dG3YsEEJCQlKSEjQ8uXLy7lyAADgjbwy8Lz++usaOXKkEhIS1KpVK02fPl1BQUGaNWtWif179uypgQMHqmXLlrrhhhs0evRotWvXTt9//305Vw4AALyR1wWeoqIirV+/XrGxsY42Hx8fxcbGKj09/Xf3NwxDqamp2r59u2655ZZS+xUWFio/P99pAwAA5nTNgefYsWOy2+2uqEWSlJeXJ5vNptDQUKf20NBQZWdnl7rfyZMnVa1aNfn7++vOO+/U22+/rV69epXaPyUlRTVr1nRsjRs3dtk5AAAA73JVgWfLli2aNGmSunbtqpCQENWrV08jRozQwoULVVBQ4Ooay6R69erKzMzU2rVr9fLLLyspKUmrVq0qtf+YMWN08uRJx3bgwIHyKxYAAJSrMgee7du366mnnlKzZs100003ae3atfrTn/6knJwcLVmyROHh4XrxxRdVt25d9enTR9OmTbuqgurWrSur1aqcnByn9pycHIWFhZV+Ij4+atq0qdq3b6+nnnpKgwYNUkpKSqn9AwICVKNGDacNAACYU5kDT1pamgoKCvTWW28pLy9PCxcu1IgRI1S3bl117txZEyZMUFZWlrZu3arbb79dixYtuqqC/P39FRUVpdTUVEeb3W5XamqqYmJiynwcu92uwsLCq6oBAACYi29ZO15Y6v17IiIiNHr0aI0ePfqqi0pKSlJ8fLyio6PVuXNnTZkyRQUFBY7nHzFihBo2bOgYwUlJSVF0dLRuuOEGFRYWasmSJfrwww+vepQJAACYS5kDT0lSUlI0ZswYZWRkqHXr1goICHBJUUOGDNGRI0c0duxYZWdnq3379lq2bJljIvP+/fvl4/O/wamCggI9+uij+uWXX1SlShW1aNFCH330kYYMGeKSegAAQMV2TYGnZ8+ekqRJkyZp8+bN8vHxUevWrdWuXTu1a9dOnTp1umS1VVklJiYqMTGxxMcunoz80ksv6aWXXrqq5wEAAOZ3TYHnwpyaTz75RJJ09uxZbdq0SRs2bNDKlSuVnJysO+64QxMmTLj2SgEAAK7SFQWehg0bKioqSlFRUerYsaM6duyohg0bOh6vUqWKOnXqpE6dOjnaoqKiCDwAAMCjrijw/L//9/+UkZGhRYsW6eWXX5bNZlNISIg6duzoFILCw8Md+/z4448uLxoAAOBKXFHgeeyxxxz/XVhYqMzMTGVkZCgjI0NLlizRa6+9puLiYp0/f97Rz8/Pz3XVAgAAXIWrnsMTEBCgLl26qGPHjlq+fLmKi4u1d+9e+fv7u7I+AACAa3ZVXy1x7tw5LV68WMOHD1dISIgSEhJktVr14Ycf6siRI66uEQAA4Jpc0QjPggULtHDhQi1dulTVq1fXwIEDtXDhQvXs2VNWq9VdNQIAAFyTKwo8w4YNU4MGDfTqq6/q4Ycflq/vNa1qB3CR8za7fK1XNfAKALiMK3pn7d69u06dOqVHH31UNWvWVExMjEaNGqVZs2YpMzPTabIygCvz456j6jBhpVbvOerpUgDAdK5oiObbb7+VJO3cuVPr1693rNCaN2+eTpw4oYCAALVt21Zr1qxxS7GAWRXb7BqzaKNOnTuvMYs3avkTt8iPkR4AcJmruibVrFkzNWvWTEOHDnW07d27V+vWrdNPP/3ksuKAyuKDtH3am1cgSdpzpEAfpO3Tw92beLgqADCPMv8JmZ2drcLCwlIfv/766zV48GBNnDhRkrRnz55rrw6oBHLyz+m1Fdud2iav2KGc/HMeqggAzKfMgeezzz5TnTp1NHDgQM2ePbvE5eerV6/Wc889p9atWysyMtKlhQJm9dJXW1RsM5zaimx2vfyvrR6qCADMp8yBJzExUVlZWerevbvmzJmjRo0a6eabb9bEiRM1cuRI1a9fXwMGDFBubq4mTZrE/XiAMkjbnacvNxyWze4ceGx2Q19kHVL6biYwA4ArXNEcnqZNmyopKUlJSUk6evSovvrqKy1ZskQRERFauHChYmJiZLFY3FUrSsAy5optUcZB+Viki/KOJMnHIi3M+EUxNwSXf2EAYDJXfSOd4OBgxcfHKz4+3pX14Ar8uOeoRv5jnd4bEa0uTfhQrIju6dhIn63/pcTH7IY0KKpROVcEAObE0EAFdfEy5mKb3dMl4SrE3BCsfu3qy+rjPDJq9bGof2QD3USQBQCXIPBUUCUtY0bF9ELfVvKzOgcef6uPnr+zpYcqAgDzIfBUQCxjNpfQGoF6uveNTm1P9W6u0BqBHqoIAMyHwFMBsYzZfOK7RqhJSFVJUpOQqorvGuHZgnDNznOZGfAqBJ4KhmXM5uRn9dHEgW1VPdBXKQPb8rUSFRzfiwZ4H95VK5gLy5hLcmEZMyqmm5oE66e/9mLFXQXHggJzYsSu4iPwVDD3dGxU4j1bJJYxmwH3VKr4WFBgPozYmQPvrhUMy5gB78WCAvNhxM48CDwVEMuYAe/EggLzYcTOPAg8FRDLmAHvw4IC82HEzlwIPBUUy5gB78KCAvNhxM5cCDwVFMuYAe/CggJzYcTOfPiUrMBYxgx4DxYUmAsjduZD4KngWMYMeA8WFJgHI3bmw6clALgICwrMgxE78yHwAIALsaDAPBixMxcCDwC4EAsKzIMRO3PhNxEAXIwFBebBiJ15EHgAwA1YUGAOjNiZh6+nCwAAwJtdGLEjxFZs/OsBAPA7CDsVH/+CAADA9Ag8AADA9Ag8AADA9Ji0XA7Gjx/v6RKuSnJysqdLANyqov5uSvx+AleKwANcoYr6IckHJMyO301cDpe0AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6Xlt4Jk6daoiIiIUGBioLl26aM2aNaX2nTlzprp3767atWurdu3aio2NvWx/AABQuXhl4FmwYIGSkpKUnJysjIwMRUZGKi4uTrm5uSX2X7VqlYYNG6ZvvvlG6enpaty4sXr37q2DBw+Wc+UAAMAbeWXgef311zVy5EglJCSoVatWmj59uoKCgjRr1qwS+3/88cd69NFH1b59e7Vo0ULvvfee7Ha7UlNTy7lyAADgjbwu8BQVFWn9+vWKjY11tPn4+Cg2Nlbp6ellOsaZM2dUXFysOnXqlNqnsLBQ+fn5ThsAADAnrws8eXl5stlsCg0NdWoPDQ1VdnZ2mY7x7LPPqkGDBk6h6WIpKSmqWbOmY2vcuPE11Q0AALyX1wWeazVp0iTNnz9fixcvVmBgYKn9xowZo5MnTzq2AwcOlGOVAACgPPl6uoCL1a1bV1arVTk5OU7tOTk5CgsLu+y+r732miZNmqSvv/5a7dq1u2zfgIAABQQEXHO9AADA+3ndCI+/v7+ioqKcJhxfmIAcExNT6n6vvPKKJkyYoGXLlik6Oro8SgUAABWE143wSFJSUpLi4+MVHR2tzp07a8qUKSooKFBCQoIkacSIEWrYsKFSUlIkSX/72980duxYzZ07VxEREY65PtWqVVO1atU8dh4AAMA7eGXgGTJkiI4cOaKxY8cqOztb7du317JlyxwTmffv3y8fn/8NTk2bNk1FRUUaNGiQ03GSk5M1bty48iwdAAB4Ia8MPJKUmJioxMTEEh9btWqV08/79u1zf0EAAKDC8ro5PAAAAK5G4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AEAAKZH4AFQ7s7b7J4uAUAlQ+ABUK5+3HNUHSas1Oo9Rz1dCoBKhMADoNwU2+was2ijTp07rzGLN6qYkR4A5YTAA6DcfJC2T3vzCiRJe44U6IO0fZ4tCEClQeABUC5y8s/ptRXbndomr9ihnPxzHqoIQGVC4AFQLl76aouKbYZTW5HNrpf/tdVDFQGoTAg8ANwubXeevtxwWDa7c+Cx2Q19kXVI6buZwAzAvQg8ANxuUcZB+VhKfszHIi3M+KV8CwJQ6RB4ALjdPR0b6aLBHQe7IQ2KalS+BQGodAg8ANwu5oZg9WtXX9aLhnmsPhb1j2ygm5oEe6gyAJWF1waeqVOnKiIiQoGBgerSpYvWrFlTat/NmzfrnnvuUUREhCwWi6ZMmVJ+hQIokxf6tpKf1Tnw+Ft99PydLT1UEYDKxCsDz4IFC5SUlKTk5GRlZGQoMjJScXFxys3NLbH/mTNn1KRJE02aNElhYWHlXC2AsgitEaine9/o1PZU7+YKrRHooYoAVCZeGXhef/11jRw5UgkJCWrVqpWmT5+uoKAgzZo1q8T+nTp10quvvqqhQ4cqICCgnKsFUFbxXSPUJKSqJKlJSFXFd43wbEEAKg2vCzxFRUVav369YmNjHW0+Pj6KjY1Venq6y56nsLBQ+fn5ThsA9/Kz+mjiwLaqHuirlIFt5Wf1urcgACblde82eXl5stlsCg0NdWoPDQ1Vdna2y54nJSVFNWvWdGyNGzd22bEBlO6mJsH66a+91IWJygDKkdcFnvIyZswYnTx50rEdOHDA0yUBlYYvIzsAypmvpwu4WN26dWW1WpWTk+PUnpOT49IJyQEBAcz3AQCgkvC6P7P8/f0VFRWl1NRUR5vdbldqaqpiYmI8WBkAAKiovG6ER5KSkpIUHx+v6Ohode7cWVOmTFFBQYESEhIkSSNGjFDDhg2VkpIi6deJzlu2bHH898GDB5WZmalq1aqpadOmHjsPAADgHbwy8AwZMkRHjhzR2LFjlZ2drfbt22vZsmWOicz79++Xj8//BqcOHTqkDh06OH5+7bXX9Nprr6lHjx5atWpVeZcPAAC8jFcGHklKTExUYmJiiY9dHGIiIiJkGKV8UQ8AAKj0vG4ODwAAgKsReAAAgOkReAAAgOkReAAAgOkReAAAgOkReAAAgOkReAAAgOkReAAAgOkReAAAgOkReAAAgOkReAAAgOkReAAAgOkReAAAgCTpvM3u6RLchsADAAD0456j6jBhpVbvOerpUtyCwAMAQCVXbLNrzKKNOnXuvMYs3qhiE470EHgAAKjkPkjbp715BZKkPUcK9EHaPs8W5AYEHgAAKrGc/HN6bcV2p7bJK3YoJ/+chypyDwIPAACV2EtfbVGxzXBqK7LZ9fK/tnqoIvcg8AAAUEml7c7TlxsOy2Z3Djw2u6Evsg4pfbd5JjATeAAAqKQWZRyUj6Xkx3ws0sKMX8q3IDci8AAAUEnd07GRLhrccbAb0qCoRuVbkBsReAAAqKRibghWv3b1Zb1omMfqY1H/yAa6qUmwhypzPQIPAACV2At9W8nP6hx4/K0+ev7Olh6qyD0IPAAAVGKhNQL1dO8bndqe6t1coTUCPVSRexB4AACo5OK7RqhJSFVJUpOQqorvGuHZgtyAwAMAQCXnZ/XRxIFtVT3QVykD28rPar544OvpAgAAgOfd1CRYP/21l3xNGHYkRngAAMB/mTXsSAQeAABQCRB4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6RF4AACA6fl6ugAA3mn8+PGeLuGqJCcne7oEwK343bw6jPAAAADTI/AAAADTI/AAAADTI/AAAADTI/AAAADTI/AAAADTI/AAAADTI/AAAADTI/AAAADTI/AAAADT89rAM3XqVEVERCgwMFBdunTRmjVrLtv/008/VYsWLRQYGKi2bdtqyZIl5VQpAADwdl4ZeBYsWKCkpCQlJycrIyNDkZGRiouLU25ubon909LSNGzYMD300EP66aefNGDAAA0YMECbNm0q58oBAIA38srA8/rrr2vkyJFKSEhQq1atNH36dAUFBWnWrFkl9n/zzTd1++2365lnnlHLli01YcIEdezYUe+88045Vw4AALyR1wWeoqIirV+/XrGxsY42Hx8fxcbGKj09vcR90tPTnfpLUlxcXKn9AQBA5eLr6QIulpeXJ5vNptDQUKf20NBQbdu2rcR9srOzS+yfnZ1d6vMUFhaqsLDQ8fPJkyclSfn5+VdbeqnOnTvn8mOWhyt9LThP78Z5XqqinqNUOc6T/2dLVlnO80qPaxjG5TsaXubgwYOGJCMtLc2p/ZlnnjE6d+5c4j5+fn7G3LlzndqmTp1q1KtXr9TnSU5ONiSxsbGxsbGxmWA7cODAZfOF143w1K1bV1arVTk5OU7tOTk5CgsLK3GfsLCwK+ovSWPGjFFSUpLjZ7vdrmPHjik4OFgWi+UazqD85Ofnq3Hjxjpw4IBq1Kjh6XLchvM0l8pwnpXhHCXO02wq6nkahqFTp06pQYMGl+3ndYHH399fUVFRSk1N1YABAyT9GkZSU1OVmJhY4j4xMTFKTU3VE0884WhbuXKlYmJiSn2egIAABQQEOLXVqlXrWsv3iBo1alSo/zmvFudpLpXhPCvDOUqcp9lUxPOsWbPm7/bxusAjSUlJSYqPj1d0dLQ6d+6sKVOmqKCgQAkJCZKkESNGqGHDhkpJSZEkjR49Wj169NDkyZN15513av78+Vq3bp1mzJjhydMAAABewisDz5AhQ3TkyBGNHTtW2dnZat++vZYtW+aYmLx//375+PxvgVnXrl01d+5cvfDCC3ruuefUrFkzff7552rTpo2nTgEAAHgRrww8kpSYmFjqJaxVq1Zd0jZ48GANHjzYzVV5l4CAACUnJ19yac5sOE9zqQznWRnOUeI8zcbs52kxjN9bxwUAAFCxed2NBwEAAFyNwAMAAEyPwAMAAEyPwAMAAEyPwFOBTZ06VREREQoMDFSXLl20Zs0aT5fkUt9995369eunBg0ayGKx6PPPP/d0SS6XkpKiTp06qXr16qpXr54GDBig7du3e7osl5s2bZratWvnuKFZTEyMli5d6umy3G7SpEmyWCxON0U1g3HjxslisThtLVq08HRZbnHw4EHdd999Cg4OVpUqVdS2bVutW7fO02W5VERExCX/nhaLRaNGjfJ0aS5F4KmgFixYoKSkJCUnJysjI0ORkZGKi4tTbm6up0tzmYKCAkVGRmrq1KmeLsVtvv32W40aNUo//vijVq5cqeLiYvXu3VsFBQWeLs2lGjVqpEmTJmn9+vVat26d/vCHP+iuu+7S5s2bPV2a26xdu1bvvvuu2rVr5+lS3KJ169Y6fPiwY/v+++89XZLLHT9+XN26dZOfn5+WLl2qLVu2aPLkyapdu7anS3OptWvXOv1brly5UpLMd6uXsnyhJ7xP586djVGjRjl+ttlsRoMGDYyUlBQPVuU+kozFixd7ugy3y83NNSQZ3377radLcbvatWsb7733nqfLcItTp04ZzZo1M1auXGn06NHDGD16tKdLcqnk5GQjMjLS02W43bPPPmvcfPPNni6j3I0ePdq44YYbDLvd7ulSXIoRngqoqKhI69evV2xsrKPNx8dHsbGxSk9P92BluFYnT56UJNWpU8fDlbiPzWbT/PnzVVBQcNnvu6vIRo0apTvvvNPpd9Rsdu7cqQYNGqhJkyYaPny49u/f7+mSXO6LL75QdHS0Bg8erHr16qlDhw6aOXOmp8tyq6KiIn300Ud68MEHK8wXaZcVgacCysvLk81mc3zVxgWhoaHKzs72UFW4Vna7XU888YS6detmyq9F2bhxo6pVq6aAgAD96U9/0uLFi9WqVStPl+Vy8+fPV0ZGhuO7/syoS5cumjNnjpYtW6Zp06Zp79696t69u06dOuXp0lxqz549mjZtmpo1a6bly5frz3/+sx5//HF98MEHni7NbT7//HOdOHFCDzzwgKdLcTmv/WoJoLIZNWqUNm3aZMq5EJJ04403KjMzUydPntRnn32m+Ph4ffvtt6YKPQcOHNDo0aO1cuVKBQYGeroct+nTp4/jv9u1a6cuXbooPDxcn3zyiR566CEPVuZadrtd0dHRmjhxoiSpQ4cO2rRpk6ZPn674+HgPV+ce77//vvr06aMGDRp4uhSXY4SnAqpbt66sVqtycnKc2nNychQWFuahqnAtEhMT9dVXX+mbb75Ro0aNPF2OW/j7+6tp06aKiopSSkqKIiMj9eabb3q6LJdav369cnNz1bFjR/n6+srX11fffvut3nrrLfn6+spms3m6RLeoVauWmjdvrl27dnm6FJeqX7/+JYG8ZcuWprx8J0k///yzvv76az388MOeLsUtCDwVkL+/v6KiopSamupos9vtSk1NNe2cCLMyDEOJiYlavHix/v3vf+v666/3dEnlxm63q7Cw0NNluNRtt92mjRs3KjMz07FFR0dr+PDhyszMlNVq9XSJbnH69Gnt3r1b9evX93QpLtWtW7dLbhOxY8cOhYeHe6gi95o9e7bq1aunO++809OluAWXtCqopKQkxcfHKzo6Wp07d9aUKVNUUFCghIQET5fmMqdPn3b6i3Hv3r3KzMxUnTp1dN1113mwMtcZNWqU5s6dq3/+85+qXr26Yw5WzZo1VaVKFQ9X5zpjxoxRnz59dN111+nUqVOaO3euVq1apeXLl3u6NJeqXr36JfOvqlatquDgYFPNy3r66afVr18/hYeH69ChQ0pOTpbVatWwYcM8XZpLPfnkk+ratasmTpyoe++9V2vWrNGMGTM0Y8YMT5fmcna7XbNnz1Z8fLx8fU0aDTy9TAxX7+233zauu+46w9/f3+jcubPx448/erokl/rmm28MSZds8fHxni7NZUo6P0nG7NmzPV2aSz344INGeHi44e/vb4SEhBi33XabsWLFCk+XVS7MuCx9yJAhRv369Q1/f3+jYcOGxpAhQ4xdu3Z5uiy3+PLLL402bdoYAQEBRosWLYwZM2Z4uiS3WL58uSHJ2L59u6dLcRuLYRiGZ6IWAABA+WAODwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CDwAAMD0CD4AKr2fPnnriiSfccuxnnnlG/fr1c8uxAZQfvloCQIV37Ngx+fn5qXr16pJ+DUDt27fXlClTrvnYJ06ckNVqdRwbQMVk0q9EBVCZ1KlTx23HrlWrltuODaD8cEkLQIXw2WefqW3btqpSpYqCg4MVGxurgoICSc6XtB544AF9++23evPNN2WxWGSxWLRv3z7Z7XalpKTo+uuvV5UqVRQZGanPPvvsss+Zl5cni8WiTZs2ufv0ALgZIzwAvN7hw4c1bNgwvfLKKxo4cKBOnTql//znPyrpivybb76pHTt2qE2bNnrxxRclSSEhIUpJSdFHH32k6dOnq1mzZvruu+903333KSQkRD169CjxebOyshQQEKAWLVq49fwAuB+BB4DXO3z4sM6fP6+7775b4eHhkqS2bduW2LdmzZry9/dXUFCQwsLCJEmFhYWaOHGivv76a8XExEiSmjRpou+//17vvvvuZQNP69at5evLWyVQ0fFbDMDrRUZG6rbbblPbtm0VFxen3r17a9CgQapdu3aZ9t+1a5fOnDmjXr16ObUXFRWpQ4cOpe6XmZmp9u3bX0vpALwEgQeA17NarVq5cqXS0tK0YsUKvf3223r++ee1evVqXX/99b+7/+nTpyVJ//rXv9SwYUOnxwICAkrdLysrSw899NC1FQ/AKzBpGUCFYLFY1K1bN40fP14//fST/P39tXjx4hL7+vv7y2azOX5u1aqVAgICtH//fjVt2tRpa9y4cYnHKCoq0tatWxUZGemW8wFQvhjhAeD1Vq9erdTUVPXu3Vv16tXT6tWrdeTIEbVs2bLE/hEREVq9erX27dunatWqqU6dOnr66af15JNPym636+abb9bJkyf1ww8/qEaNGoqPj7/kGFu3blVxcTGBBzAJAg8Ar1ejRg199913mjJlivLz8xUeHq7JkyerT58+JfZ/+umnFR8fr1atWuns2bPau3evJkyY4FittWfPHtWqVUsdO3bUc889V+IxMjMzFR4ezn14AJPgTssAUILExETl5ubqk08+8XQpAFyAOTwA8Bvnzp3T+vXrtXDhQsXFxXm6HAAuQuABgN+YMmWKYmNjddddd2nEiBGeLgeAi3BJCwAAmB4jPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPQIPAAAwPT+P6sCxDcqrFYMAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(range(N_s), occs[0], color='grey', label='Exact')\n",
    "plt.scatter(range(N_s), qc_occ_nums[0], marker='d', label='QC')\n",
    "plt.legend()\n",
    "plt.xlabel(r'site $j$')\n",
    "plt.ylabel(r'$\\langle N_j \\rangle$')\n",
    "title = r'Initial State: '\n",
    "title += r'$N_s = $' + str(N_s)\n",
    "title += r', $n_k = \\pm$' + str(nk_mean)\n",
    "title += r', $\\sigma_k=$' + str(sig_k)\n",
    "plt.title(title)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "68.57814034807876"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "perc_errors = []\n",
    "for i in range(N_s):\n",
    "    err = np.abs(1 - qc_occ_nums[0][i]/occs[0][i]) * 100\n",
    "    perc_errors.append(err)\n",
    "np.mean(perc_errors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[12.868061552895615,\n",
       " 2.4812273111804894,\n",
       " 17.894319712510253,\n",
       " 352.9533083400625,\n",
       " 9.582965698137102,\n",
       " 2.125269415643616,\n",
       " 8.700801101912049,\n",
       " 142.01916965228838]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "perc_errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f233f9996d0>"
      ]
     },
     "execution_count": 152,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhIAAAHNCAYAAABLvZLYAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQipJREFUeJzt3XlcVnX6//H3DcpiCC4kqKG4VGoqmqSp0yoT46TWd0bTxpKobJPUYWzMX6Vlk9jmYOVXs0lzWu3bjJPVaBm55ORSbpmNmqVJGqCVgCig3Of3h8M93QJ63+fch8MNr+c8zmPiw1mucxC4uD7LcRmGYQgAAMCEEKcDAAAAwYtEAgAAmEYiAQAATCORAAAAppFIAAAA00gkAACAaSQSAADANBIJAABgGokEAAAwjUQCAACYRiJRhz3xxBPq0qWL3G63X8fNmzdP7dq1U1lZmU2R1X3+PDueFwCYRyJRRxUVFenxxx/X5MmTFRLi35fplltuUXl5uZ5//nmf9n/ppZfkcrlq3NavX1/tfhEREWrTpo1SU1P1zDPPqLi42O/7XLRokZo0aaJjx4552t5//33PNb788ssqxwwdOlTnnXdejec807MzDEPTp0/Xxx9/7Gk70/N64IEH1LFjR7/v63RW78lOX331lUaNGqXzzjtPTZo0UZcuXTR9+nSvr4kZb7/9tlwul/7v//4vQJFWVVZWpsmTJ6tNmzaKjIxUv379tGLFClPneuyxx+RyudS9e/cAR1l7rDyPHTt2aMSIEerYsaOaNGmi2NhYXX755XrnnXdsjhrBjkSijlqwYIFOnjypG2+80e9jIyIilJaWplmzZsmfd7JNnz5dL7/8cpWtc+fO1e43d+5c3XvvvZKkiRMnqkePHvr888/9inXp0qUaNGiQmjRp4mnbtm2bJCkkJETvvvtulWO2bdumnj171njOMz273bt3a9q0afr+++89bWd6XkOHDtXevXu1Y8cOv+6rupgl8/dkl9zcXPXt21fr169XRkaGsrOz1b9/f02bNs3Uv72f2759uyTZel+33HKLZs2apdGjR2v27NkKDQ3Vr3/9a61du9av83z33XeaMWOGzjnnHJsirR1Wnse3336r4uJipaWlafbs2XrooYckScOGDdP8+fPtDh3BzECd1LNnT+Omm24yffxnn31mSDJycnLOuu/ChQsNScann35qer+cnBwjMjLSaN++vXHs2DGfYiwtLTWioqKM+fPne7WPHj3aiI6ONgYPHmz84he/8Prcjz/+aEgy7r///hrPe6Zn9+qrrxqSjK+++sqrvabn5Xa7jbi4OCMrK8une6qJ1Xuyy2OPPWZIMr744guv9jFjxhiSjB9//NH0uW+44QYjMjLSOHnypNUwq7VhwwZDkvHkk0962o4fP2506tTJ6N+/v1/nGjlypHH11VcbV1xxhXHRRRcFOtRaEcjnUenkyZNGUlKSceGFFwYqTNRDVCTqoL179+rzzz9XSkqKV/u3336re+65RxdeeKEiIyPVsmVLjRgxQvv27atyjj59+qhFixZ6++23ayXmq6++Wg899JC+/fZbvfLKKz4ds3LlSpWUlGjo0KFe7du2bVOPHj00ZMgQrVu3Tj/88IPX56Sa/8qt6dlJUt++fTV69GhJ0vnnny+Xy6VmzZpJqvl5uVwuDRkypNry7s6dO7V//36f7tXKPdmpqKhIkhQXF+fV3rp1a4WEhCgsLMz0ubdv365u3bopNDTUUow1eeuttxQaGqo77rjD0xYREaHbbrtN69atU25urk/nWbNmjd566y1lZ2cHLLZFixapZ8+eioyMVJ8+fbRu3TrdcMMNSkpKCtg1Theo5/FzoaGhSkhI0JEjRwIYKeobEok66JNPPpEkXXzxxV7tn376qT755BONGjVKzzzzjO666y7l5OToyiuvrLY/++KLL9a//vUvn69bWFiow4cPe20//4V3NjfffLMk6YMPPvBp/3feeUeXXHKJ4uPjPW3l5eXatWuXkpKSNGTIEFVUVOif//yn5/OVXSc1/UCu6dlJ0uTJk9W9e3edf/75nm6bnyc9NT2vYcOGaf369Tp06JBXe9euXTVmzJiz3qfVe6rJiRMnqny9atpqGnR65ZVXSpJuu+02bd26Vbm5uVq8eLHmzp2r8ePHmy71l5WV6auvvlKPHj1si33Lli264IILFB0d7dXet29fSdLWrVvPGmdFRYXuvfde3X777TXG6q+srCzdcsst6ty5s5555hldeumlGjZsmDZv3lznn4cklZSU6PDhw/r666/15z//WcuWLdOgQYN8fwBoeJwuiaCqBx980JBkFBcXe7VX12Wwbt06Q5Lx17/+tcrn7rjjDiMyMvKs16vssqhuCw8Pr7LfmbpAYmJijN69e5/1moZhGAkJCcajjz7q1bZlyxZDkjFv3jzDMAyjR48exogRIzyfv/XWW43w8PAay+U1PbtK7dq1M2655ZZqP1fT8yopKTEiIyONl156yatdknHFFVfUeH+BuqearFy5ssav2+nb3r17azzPo48+akRGRnrt/8ADD/gVy+kq73nWrFm2xX7RRRcZV199dZX2HTt2eD3vM3nuueeMmJgYo6CgwDAMw3LXxqZNm4zQ0FBjypQpXu233XabIanGLrK68jwMwzDuvPNOz7VCQkKM4cOHW+riQv3XKLBpCQLhhx9+UKNGjRQVFeXVHhkZ6fnvEydOqKioSJ07d1azZs20efNmT0WgUvPmzXX8+HEdO3bMazBjTebMmaMLLrjAq83fsnRUVJRPsze2bNmi3NxcDRs2zKu98q/zyjL/kCFDNGfOHJ04cUKNGzfWtm3bdNFFF3nicrvdSkhI0ObNmxUXF1fjs5NOVVz2799fYxdCTc+rSZMmGjRokN555x2lpaV52g0fB7L6e0++SkpK8nlE/s+rPqdLTEzU5Zdfrt/+9rdq2bKl3nvvPc2YMUPx8fHKyMjwK6ZKlQMta/oLPBCxHz9+XOHh4VXaIyIiPJ8/kx9++EFTp07VQw89pHPPPdenWM5mxowZat26tR555BGv9ssuu0wvvvhinX4elSZOnKjhw4fr4MGDevPNN1VRUaHy8nKfjkXDRCIRRI4fP66srCwtXLhQBw4c8PpFVlhYWGX/ys+7XC6fzt+3b18lJydbivHo0aNq1arVWfdbvXq1mjVrVuWX+rZt2+RyuTw/cIcMGaKsrCytWbNGV155pXbs2KFRo0Z59g8JCdGBAwd8iu30X+inO9PzuvLKK5WVleXTdU7n7z35qnnz5tWOBfHHG2+8oTvuuEO7d+/2TD/9zW9+I7fbrcmTJ+vGG29Uy5Yt/T7v2WZsBCL2yMjIatf+KC0t9Xz+TB588EG1aNHCM/PIqrKyMi1btkx33323Gjdu7PW5o0ePSqo5saoLz6NSly5d1KVLF0nSmDFjdM0112jo0KHasGGDzz9L0LCQSNRBLVu21MmTJ1VcXKymTZt62u+9914tXLhQEydOVP/+/RUTEyOXy6VRo0ZV22/6008/qUmTJj7/ALHqu+++U2FhYZXpotXp06ePjhw5oj179njt//nnn6tjx46eisKll16q2NhYvfPOO2rTpo1KS0vPOJagpmdXeW6p5rEIZ3peGzduNJ1kWb2nmpSXl+vHH3/0ad9zzz232orH//7v/6p3795V1rAYNmyYXnrpJW3ZssXUL7jt27erVatWNSaVgYi9devW1SaRlVN727RpU+M5v/rqK82fP1/Z2dk6ePCgp720tFQnTpzQvn37FB0drRYtWvgUoyR9/fXXOnbsWLVfy2+++UbR0dFq165dtcc6/TzOZPjw4brzzju1e/duXXjhhabOgfqNRKIOqvxrYO/evV5/0b311ltKS0vT008/7WkrLS2tcUT13r171bVrV1tj/bmXX35ZkpSamnrWfQcMGKDY2FgtXbpUmZmZnvbPP/9cAwcO9HwcEhKiwYMH65133tGll14qyfuv3Pnz52vVqlV67bXXJNX87CrP3bp1a8XGxlYbU03P68SJE1q+fLlmzJhx1vuqjr/3tHPnTmVkZGjz5s0yDEM333yznnnmmSrn/eSTT3TVVVf5FMPevXuVmJhYpT0/P1/Nmzev0n7ixAlJ0smTJ306/+m2b99+xlkogYi9V69eWrlypYqKirwGGG7YsMHz+ZocOHBAbrdb48eP1/jx46t8vkOHDpowYYJfMzkquw5OXwStrKxMr7/++hkXunL6eZxJ5X1VV/UEJBKJOql///6SpM8++8zrh3FoaGiVfvlnn31WFRUV1Z5n8+bNnumOdvvoo4/06KOPqkOHDj5ds3KhnHfeeceTSOTl5amgoKDKX3RDhgzRyy+/rNdff12Sd0Vhx44dXj+ga3p2krR///4zrh5Z0/NavXq1ioqKqkxT3blzp5o0aVLjX5lm72n06NGaPHmyRowYoeLiYn311VfVnjsQ/eoXXHCBPvjgA+3evdtrfMzrr7+ukJAQU1NSf/rpJx04cEA33HBDjfsEIvbhw4frqaee0vz58zVp0iRJp35pL1y4UP369VNCQoIk6dixY9q/f79iY2M9SWT37t21ZMmSKud88MEHVVxcrNmzZ6tTp04+xVepffv2kqS1a9d6/TuaOXOmvv/++yr/fn7O6echSQUFBVUqSCdOnNBf//pXRUZGqlu3bj7Fh4aHRKIO6tixo7p3764PP/xQt956q6e98pdPTEyMunXrpnXr1unDDz+stg9706ZN+vHHH3Xdddf5fN1ly5Zp586dVdoHDBjgtUx05X4nT55Ufn6+PvroI61YsULt27fX0qVLPYO7zmbYsGEaNWqUfvrpJzVv3rzG9RRSU1PVuHFjT1fAz+93x44duvrqqz0f1/TspFN/ZX700Ud64okn1KZNG3Xt2lV9+vSRdObntXTpUiUlJVVJGLp27aorrrhCq1atqvEezdzT119/rfLycrndbkVHR3tiPF0g+tXvu+8+LVu2TJdddpkyMjLUsmVLvfvuu1q2bJluv/32KuVwl8t11nv2ZUXLQMTer18/jRgxQlOmTFFBQYE6d+6sRYsWad++fXrxxRc9+23cuFFXXXWVpk2bpocffliSFBsbq+uvv77KOSsrEKd/zpf7jo2N1a9+9SvNnz9fjRo1Us+ePfXhhx9q06ZNkmoeHyE5/zwk6c4771RRUZEuv/xytW3bVnl5eXr11Ve1c+dOPf3009UOYAYkMf2zrpo1a5YRFRXlNeXzp59+MtLT043Y2FgjKirKSE1NNXbu3Gm0b9/eSEtL8zp+8uTJRrt27Qy3233Wa51p+qckY+HChdXuFxYWZsTHxxu//OUvjdmzZxtFRUV+3WNRUZERFhZmvPrqq4ZhGMYTTzxhSDL27NlTZd+rrrrKkGT86le/8mqPj4+vskpldc/OMAzjwIEDRmpqqhEVFWVIMp555hnP5870vBITE40HH3ywSrt8mP5p5p7++c9/GgMHDjTi4uKM++67zzhx4sQZr2HVhg0bjMGDBxvx8fFG48aNjQsuuMB47LHHqly3uLjYkGSMGjXqjOd77rnnDEnG5s2b7QzbMIxTKzdOmjTJiI+PN8LDw41LLrnEWL58udc+lVMrp02bdtbzVTf909f7NgzDyM/PN66//nojOjraaN26tTFhwgRjyZIlhiRj/fr1ft2bGVaex+uvv26kpKQYcXFxRqNGjYzmzZsbKSkpxttvv2173AhuJBJ11JEjR4wWLVoYf/nLX/w+trS01IiPjzeys7NtiCywUlNTjZEjR5o69ocffjAiIyONiooKr3Z/n92Zntfnn39uSDI2bNhgKkYr9u3bZ7Rr1854//33a/3a1XnvvfcMl8tlfP7552fc73e/+53RtGlT2xOg2uLrfddk3LhxRmxsbL15HsDpWNmyjoqJidEf//hHPfnkk36/RnzhwoVq3Lix7rrrLpuiC5xhw4Zp+fLlnsF9/tixY4e6du1aZXCbv8/uTM9r6dKlio+P1yWXXOJ3fGb87W9/0969eyWdGmtQXl7uGUDqtJUrV2rUqFFnLNEfOnRIS5cu1fDhw9WoUf3oOfXlvqVT4xGM08YwrVmzRs8//7zuuuuuevM8gNO5jNP/5QO1qLS0VN99950SExP9/kE7b948ffLJJ/rrX/9qU3SnZjVUVFSYnjrnr/Hjx+vNN9/U0aNH1alTJ82YMUPXXnttrVzbiu3btysnJ0dz585VXl6etm7dqg4dOjgdVq1aunSppk+fruHDh3sWiVu4cKF69+6t1atX19o0bKC2kUggaN1zzz3q1KmT/vCHPzgdSoOXlpamv//97xo4cKAef/xxW19OVVetXbtWf/zjH/Xll1/q2LFjSkxM1A033KApU6YE/evJgTMhkUBQOnr0qHr06KE333yz1rodAABVMUYCQWfDhg264IIL9Jvf/IYkAgAcRkUCAACYRkUCAACYFvTzkdxutw4ePKimTZvyZjoAwBkZhqHi4mK1adOmytTxQCotLQ3I69fDwsJ8Xi3YKUGfSBw8eNCzhjwAAL7Izc0947t3rCgtLVWHDvHKy7P+orP4+Hjt3bu3TicTQZ9IVL4qes01fRTVOLhv59Ovz3c6BMs+Lqgf09y+O+7/All10XmRjZ0OwbLLWpU4HUJAXNKp+pevBZN2E4udDsGy4pIKdbh+l+d3hx3Ky8uVl1eofbmzFR1tfv2QoqLjSkyYoPLychIJO1V2Z0Q1bqSmQZ5IRIaGOR2CZWEh4U6HEBCN6kk3WVhI8P+bigytH0ldsP98kqToc0KdDiFgaqMrPCoqXFFR5n8m+ruqsVOC/182AAB1kGGclGGctHR8MCCRAADABoZRIcOosHR8MGD6JwAAMI2KBAAANnAbJ+W20D1h5djaRCIBAIANGsoYCbo2AACAaVQkAACwwanBllYqEsEx2JJEAgAAGxjukzLcFhIJC8fWJro2AACAaVQkAACwg3Hy1Gbl+CBAIgEAgA2YtQEAAHAWVCQAALCD+6TktvDSuSAZbEkiAQCADU51bZh/Y2qwdG2QSAAAYAf3Sclt4dXrQVKRYIwEAAAwjYoEAAB2aCAVCRIJAABsUWFxLYjgWCKbrg0AAGAaFQkAAGzgcp+Uy23+73UXXRsAADRg7pOShUQiWMZI0LUBAABMc7wiUVFRoYcfflivvPKK8vLy1KZNG91yyy168MEH5XK5fD5P48Yn1bixYWOk9vuhLNzpECw7cDw4MuizSYgMczqEgKgPX4/68H0hSY0bW1jhsI4o7X690yFYVlpUJunL2rlYA6lIOJ5IPP7445o7d64WLVqkiy66SJ999pnS09MVExOj8ePHOx0eAACmuIyTchkWxkiwsqVvPvnkE1133XW69tprJUmJiYl6/fXXtXHjRocjAwAAZ+P4GIkBAwYoJydHu3fvliRt27ZNa9eu1eDBg6vdv6ysTEVFRV4bAAB1jtstuSssbG6n78Anjlck7r//fhUVFalLly4KDQ1VRUWFHnvsMY0ePbra/bOysvTII4/UcpQAAPjn1PRP38f6VXd8MHC8IvHmm2/q1Vdf1WuvvabNmzdr0aJFeuqpp7Ro0aJq958yZYoKCws9W25ubi1HDACADyxVI/6zBQHHKxL33Xef7r//fo0aNUqS1KNHD3377bfKyspSWlpalf3Dw8MVHl4/RnEDABDsHE8kjh07ppAQ78JIaGio3EHSNwQAQLXcJyULXRtM//TR0KFD9dhjj6ldu3a66KKLtGXLFs2aNUu33nqr06EBAGCay11hcYlsujZ88uyzz+qhhx7SPffco4KCArVp00Z33nmnpk6d6nRoAADgLBxPJJo2bars7GxlZ2c7HQoAAIFjVFhb2dKgIgEAQIPlcrstdU+4gmSsoOPTPwEAQPCiIgEAgB3cFRZnbdC1AQBAg3Vq1oaVlS2DI5GgawMAAJhGRQIAADs0kK4NKhIAANjgVNeGtc2MOXPmKDExUREREerXr582btxY474vvfSSXC6X1xYREeHX9ahIAABgBwcqEosXL1ZmZqbmzZunfv36KTs7W6mpqdq1a5datWpV7THR0dHatWuX52OXy7+YqUgAAFBPzJo1S2PHjlV6erq6deumefPmqUmTJlqwYEGNx7hcLsXHx3u2uLg4v65JIgEAgA1cbuM/i1KZ3QxJUlFRkddWVlZW7fXKy8u1adMmpaSkeNpCQkKUkpKidevW1Rjn0aNH1b59eyUkJOi6667Tjh07/LpPEgkAAOzgrrC+SUpISFBMTIxny8rKqvZyhw8fVkVFRZWKQlxcnPLy8qo95sILL9SCBQv09ttv65VXXpHb7daAAQP03Xff+XybjJEAAKAOy83NVXR0tOfj8PDwgJ27f//+6t+/v+fjAQMGqGvXrnr++ef16KOP+nQOEgkAAOxgVEhWXpfxn5d2RUdHeyUSNYmNjVVoaKjy8/O92vPz8xUfH+/TJRs3bqzevXtrz549PodZbxKJsuPhanyy3txO0CowipwOISA6NWrhdAgBUT++Hk2dDiAgSo/7N6WuLmr87cdOh2BZ46Mna+1aLsMtl2FhZUvDvywkLCxMffr0UU5Ojq6//npJktvtVk5OjjIyMnw6R0VFhbZv365f//rXPl+X37wAANQTmZmZSktLU3Jysvr27avs7GyVlJQoPT1dkjRmzBi1bdvWM85i+vTpuvTSS9W5c2cdOXJETz75pL799lvdfvvtPl+TRAIAADu4LXZtmFhHYuTIkTp06JCmTp2qvLw89erVS8uXL/cMwNy/f79CQv47z+Knn37S2LFjlZeXp+bNm6tPnz765JNP1K1bN5+vSSIBAIAd3G6LC1KZy0IyMjJq7MpYtWqV18d//vOf9ec//9nUdSox/RMAAJhGRQIAADs4VJGobSQSAADY4NTqlNaODwYkEgAA2MHttjjYMjgSCcZIAAAA06hIAABghwZSkSCRAADADg0kkaBrAwAAmEZFAgAAOxgVktuwcHxwVCRIJAAAsEFDmf5J1wYAADCNigQAAHZoIIMtSSQAALBDA0kk6NoAAACmUZEAAMAObsNaVcHKjI9aRCIBAIAd3IbFrg0SCQAAGi7LrxEPjkSCMRIAAMA0KhIAANihgVQkSCQAALBDAxkjQdcGAAAwrd5UJI6XRSi0Irhvp0fLQ06HYFliQbzTIQREiIVqZF2S2KiZ0yFY1qNlntMhBMTx0ginQ7AsZMXHTodgWUhpLf6Vb7glw8IPEyM4KhLB/ZsXAIC6yrDYtREkiQRdGwAAwDQqEgAA2KGBDLYkkQAAwA4NJJGgawMAAJhGRQIAABsY7lObleODAYkEAAB2aCBdGyQSAADYwS2LiUSgArEXYyQAAIBpVCQAALBDA6lIkEgAAGAH4z+bleODAF0bAADANCoSAADYwHC7ZLjNv7SL6Z8AADRkDWSMBF0bAADANCoSAADYwXBJFro2gmWwJYkEAAA2aChjJOjaAAAAplGRAADADm6LXRtBUpEgkQAAwA6G69Rm+vjAhWInEgkAAGzAGAkAAICzoCIBAIAd3CEWx0gER98GiQQAAHZgsGVwyS9sruLQxk6HYcmFCfudDsGyXxw+1+kQAqLohNMRBEaHcyqcDsGyDm0OOh1CQOzKbed0CJYZ713udAiWHT15UlKO02HUK/UmkQAAoC4xDJcMC7M2jODo2SCRAADAFg1kjASzNgAAgGlUJAAAsIHhlsV1JIKjIkEiAQCAHSy//dPCsbWIrg0AAGAaFQkAAGxgfdYGFQmfHThwQDfddJNatmypyMhI9ejRQ5999pnTYQEAYJ47xPoWBByvSPz0008aOHCgrrrqKi1btkznnnuuvvrqKzVv3tzp0AAAMM36S7uCoyLheCLx+OOPKyEhQQsXLvS0dejQwcGIAACArxyvmyxdulTJyckaMWKEWrVqpd69e+uFF16ocf+ysjIVFRV5bQAA1DWVYySsbMHA8UTim2++0dy5c3X++efr/fff1913363x48dr0aJF1e6flZWlmJgYz5aQkFDLEQMA4IMGMkbC8SjdbrcuvvhizZgxQ71799Ydd9yhsWPHat68edXuP2XKFBUWFnq23NzcWo4YAIC6a86cOUpMTFRERIT69eunjRs3+nTcG2+8IZfLpeuvv96v6zmeSLRu3VrdunXzauvatav276/+TZjh4eGKjo722gAAqGsqB1ta2fy1ePFiZWZmatq0adq8ebOSkpKUmpqqgoKCMx63b98+TZo0SZdddpnf13Q8kRg4cKB27drl1bZ79261b9/eoYgAALDOiTESs2bN0tixY5Wenq5u3bpp3rx5atKkiRYsWFDjMRUVFRo9erQeeeQRdezY0e9rOp5I/P73v9f69es1Y8YM7dmzR6+99prmz5+vcePGOR0aAABBo7y8XJs2bVJKSoqnLSQkRCkpKVq3bl2Nx02fPl2tWrXSbbfdZuq6jk//vOSSS7RkyRJNmTJF06dPV4cOHZSdna3Ro0c7HRoAAOZZHTDpPvV/p89ODA8PV3h4eJXdDx8+rIqKCsXFxXm1x8XFaefOndVeYu3atXrxxRe1detW02E6nkhI0pAhQzRkyBCnwwAAIGACtSDV6bMTp02bpocffthKaJKk4uJi3XzzzXrhhRcUGxtr+jx1IpEAAADVy83N9ZpYUF01QpJiY2MVGhqq/Px8r/b8/HzFx8dX2f/rr7/Wvn37NHToUE+b232qDNKoUSPt2rVLnTp1Omt8JBIAANggUC/t8nWGYlhYmPr06aOcnBzPFE63262cnBxlZGRU2b9Lly7avn27V9uDDz6o4uJizZ492+d1mkgkAACwg2FxjITh/yGZmZlKS0tTcnKy+vbtq+zsbJWUlCg9PV2SNGbMGLVt21ZZWVmKiIhQ9+7dvY5v1qyZJFVpPxMSCQAAbODES7tGjhypQ4cOaerUqcrLy1OvXr20fPlyzwDM/fv3KyQksBM2SSQAAKhHMjIyqu3KkKRVq1ad8diXXnrJ7+vVm0Ri+48tFBES5nQYllzU8WunQ7DsvHNKnA4hIP55oJnTIQTEr9sG/0vtGjU+4XQIAbHth5ZOh2DZRwdbOR2CZWXuMkk5tXItw5DFMRIBDMZG9SaRAACgTrHYtSErx9Yix1e2BAAAwYuKBAAANjCMEBmG+b/XjSDp2yCRAADADm6Xte4JujYAAEB9R0UCAAAbBGply7qORAIAABs4sSCVE+jaAAAAplGRAADABszaAAAApjWUrg0SCQAAbNBQBlsyRgIAAJhGRQIAABs0lIoEiQQAADYwDItjJIIkkaBrAwAAmEZFAgAAGzD9EwAAmNZQpn/StQEAAEyjIgEAgA2YtQEAAEwjkQAAAKYZbmvjHAx3AIOxEWMkAACAaVQkAACwAV0bAADANOvrSARHp0G9SSR2F4UoLCTU6TAsKT0e4XQIll2WtM3pEAKioPRSp0MIiPrw9Sg+Eu10CAGxsyi4fz5J0o7jhU6HYFmFUe50CPVOvUkkAACoS9yGS24L3RNWjq1NJBIAANjB4sqWYmVLAABQ31GRAADABszaAAAApjWURIKuDQAAYBoVCQAAbNBQKhIkEgAA2MBthMhtYVEpK8fWJhIJAABsYBjWpn8GS0UiONIdAABQJ1GRAADABoyRAAAApjWURIKuDQAAYBoVCQAAbMBLuwAAgGl0bQAAAJwFFQkAAGxAReIMBgwYoKKiokDHAgBAvVE5RsLKFgxMJRLr169XaWlplfaioiJNnjzZclAAACA4+JVIDB8+XDNnzpTL5VJBQUGVz5eUlOipp54KWHAAAAQrw/hv94a5zek78I1fYyTatWund999V4ZhKCkpSS1btlRSUpKSkpLUq1cv7dq1S61bt7YrVgAAgkZDGSPhVyIxa9YsSVJYWJj+9a9/6eDBg9qyZYu2bt2qJUuWyO1264knnrAlUAAAgolhcZxDvUwkKpWUlKhx48aSpOuuuy6gAQEAgOBhKpGoTCIAAED16NoIMnmlJ9XIFdzrax08fK7TIVh22Q2fOR1CQFxzJNrpEAKi1TW7nQ7Bsl3/O9jpEALi4PETTodg2b8r/uV0CJYZhrsWr9UwEong/s0LAAAcVW8qEgAA1CW8tAsAAJhG18ZZfPzxx7rpppvUv39/HThwQJL08ssva+3atQELDgAA1G2mEom//e1vSk1NVWRkpLZs2aKysjJJUmFhoWbMmBHQAAEACEa8a+MM/vSnP2nevHl64YUXvKaCDhw4UJs3bw5YcAAABCtDLstbMDCVSOzatUuXX355lfaYmBgdOXLEakwAACBImEok4uPjtWfPnirta9euVceOHS0HBQBAsLP2wi5rAzVrk6lEYuzYsZowYYI2bNggl8ulgwcP6tVXX9WkSZN09913BzpGAACCDmMkzuD+++/X7373Ow0aNEhHjx7V5Zdfrttvv1133nmn7r333kDHCABA0HGqIjFnzhwlJiYqIiJC/fr108aNG2vc9+9//7uSk5PVrFkznXPOOerVq5defvllv65nah0Jl8ulBx54QPfdd5/27Nmjo0ePqlu3boqKijJzOgAAEACLFy9WZmam5s2bp379+ik7O1upqanatWuXWrVqVWX/Fi1a6IEHHlCXLl0UFhamd999V+np6WrVqpVSU1N9uqbpBalKS0v1+eefq6CgQG63W3l5eZ7PDRs2zOxpAQCoF9yyuLKliVkbs2bN0tixY5Weni5Jmjdvnt577z0tWLBA999/f5X9r7zySq+PJ0yYoEWLFmnt2rX2JhLLly/XzTffrB9++KHK51wulyoqKsycFgCAeiNQK1sWFRV5tYeHhys8PLzK/uXl5dq0aZOmTJniaQsJCVFKSorWrVvnw/UMffTRR9q1a5cef/xxn+M0NUbi3nvv1Q033KDvv/9ebrfbayOJAAAgcBISEhQTE+PZsrKyqt3v8OHDqqioUFxcnFd7XFycV6/B6QoLCxUVFaWwsDBde+21evbZZ/XLX/7S5/hMVSTy8/OVmZlZJVgAAHCKWy5T3RM/P16ScnNzFR0d7WmvrhphRdOmTbV161YdPXpUOTk5yszMVMeOHat0e9TEVEVi+PDhWrVqlZlDz2jmzJlyuVyaOHFiwM8NAECtsjpj4z9dG9HR0V5bTYlEbGysQkNDlZ+f79Wen5+v+Pj4GsMMCQlR586d1atXL/3hD3/Q8OHDa6x6VMdUReK5557TiBEj9PHHH6tHjx5ey2RL0vjx4/0+56effqrnn39ePXv2NBMSAAANWlhYmPr06aOcnBxdf/31kiS3262cnBxlZGT4fB632+15h5YvTCUSr7/+uj744ANFRERo1apVcrn+W7pxuVx+JxJHjx7V6NGj9cILL+hPf/qTmZAAAKhTrC4qZebYzMxMpaWlKTk5WX379lV2drZKSko8szjGjBmjtm3beioOWVlZSk5OVqdOnVRWVqZ//vOfevnllzV37lyfr2kqkXjggQf0yCOP6P7771dIiOk3kXuMGzdO1157rVJSUs6aSJSVlXllSqePZgUAoC4I1KwNf4wcOVKHDh3S1KlTlZeXp169emn58uWeMY379+/3+r1dUlKie+65R999950iIyPVpUsXvfLKKxo5cqTP1zSVSJSXl2vkyJEBSSLeeOMNbd68WZ9++qlP+2dlZemRRx6xfF0AAOqjjIyMGrsyTh/f+Kc//clyT4CpTCAtLU2LFy+2dGHp1EjUCRMm6NVXX1VERIRPx0yZMkWFhYWeLTc313IcAAAEmjsAWzAwVZGoqKjQE088offff189e/asMthy1qxZPp1n06ZNKigo0MUXX+x17jVr1ui5555TWVmZQkNDvY6paSEOAADqEie6NpxgKpHYvn27evfuLUn64osvvD7384GXZzNo0CBt377dqy09PV1dunTR5MmTqyQRAAAEC7dhbsDkz48PBqYSiZUrVwbk4k2bNlX37t292s455xy1bNmySvvZHHD9oFBX47PvWIdtO1z1hSrB5pIB1zgdQkC037/K6RACorQefD22TQ/+7wtJ+s512OkQLCs7+aPTIVhmGEHy2zmImH5pFwAAqJkhlwwLK1taObY2+ZxIZGZm6tFHH9U555yjzMzMM+7r6xiJ6tixYiYAALXNiXUknOBzIrFlyxadOHHC89818WeMBAAACG4+JxI/HxexaNEinXfeeVXWkTAMg+mYAACocrClteODgal1JDp06KDDh6sOHPrxxx/VoUMHy0EBABDsKsdIWNmCgalEoqZRr0ePHvV5YSkAABD8/Jq1UTnI0uVyaerUqWrSpInncxUVFdqwYYN69eoV0AABAAhGDLasRuUgS8MwtH37doWFhXk+FxYWpqSkJE2aNCmwEQIAEIQM49Rm5fhg4FciUTngMj09XbNnz1Z0dLQtQQEAgOBgakGqhQsXBjoOAADqFUMuuVmQCgAAmMFLuwAAgGkNZbClqemfAAAAEhUJAABsYfxns3J8MCCRAADABnRtAAAAnAUVCQAAbOD+z2bl+GBAIgEAgA0ayvRPujYAAIBpVCQAALBBQxlsSSIBAIANGsr0T7o2AACAaVQkAACwAV0bAADANKZ/BplCV4FCXMF9O3uPtnY6BMtCElKdDiEgGnX92OkQAsJdD74ee4/udDqEgDjiync6BMvcxkmnQ7DMMGpv5AHTPwEAAM4iuP+EBwCgjjJkrXsiWGZtkEgAAGADQxa7NkTXBgAAqOeoSAAAYAO3cWqzcnwwIJEAAMAGrGwJAABwFlQkAACwAStbAgAA0xrKypZ0bQAAANOoSAAAYIOGskQ2iQQAADZoKF0bJBIAANjAME5tVo4PBoyRAAAAplGRAADABm655Lbwvgwrx9YmEgkAAGzQUJbIpmsDAACYRkUCAAA7WBxsGSwv2yCRAADABg1ljARdGwAAwDQqEgAA2KChrCNBIgEAgA0aysqWdG0AAADTqEgAAGCDhrKORL1JJMrcxXK5Qp0Ow5KW4cFSyKpZWONop0MIiIp9x50OISDqw9ejPnxfSFJpYZHTIaCWGbI2gzNI8oj6k0gAAFCXnKpIWJj+GSSZBGMkAACAaVQkAACwAdM/AQCAaUz/BAAAOAsSCQAAbFDZtWFlM2POnDlKTExURESE+vXrp40bN9a47wsvvKDLLrtMzZs3V/PmzZWSknLG/atDIgEAgA3cAdj8tXjxYmVmZmratGnavHmzkpKSlJqaqoKCgmr3X7VqlW688UatXLlS69atU0JCgq655hodOHDA52uSSAAAUE/MmjVLY8eOVXp6urp166Z58+apSZMmWrBgQbX7v/rqq7rnnnvUq1cvdenSRX/5y1/kdruVk5Pj8zVJJAAAsIFh/Hd1SzNbZddGUVGR11ZWVlbt9crLy7Vp0yalpKR42kJCQpSSkqJ169b5FPOxY8d04sQJtWjRwuf7JJEAAMAGRgA2SUpISFBMTIxny8rKqvZ6hw8fVkVFheLi4rza4+LilJeX51PMkydPVps2bbySkbNh+icAAHVYbm6uoqP/u9x9eHi4LdeZOXOm3njjDa1atUoRERE+H0ciAQCADQL10q7o6GivRKImsbGxCg0NVX5+vld7fn6+4uPjz3jsU089pZkzZ+rDDz9Uz549/YqTrg0AAGxQ29M/w8LC1KdPH6+BkpUDJ/v371/jcU888YQeffRRLV++XMnJyX7fJxUJAABs4MTKlpmZmUpLS1NycrL69u2r7OxslZSUKD09XZI0ZswYtW3b1jPO4vHHH9fUqVP12muvKTEx0TOWIioqSlFRUT5dk0QCAIB6YuTIkTp06JCmTp2qvLw89erVS8uXL/cMwNy/f79CQv7bGTF37lyVl5dr+PDhXueZNm2aHn74YZ+uSSIBAIANAjVGwl8ZGRnKyMio9nOrVq3y+njfvn3mLvIzJBIAANjg51M4zR4fDBhsCQAATKMiAQCADZzq2qhtJBIAANjAyhs8K48PBnRtAAAA06hIAABgAyfWkXACiQQAADZwy+IYiYBFYi/HE4msrCz9/e9/186dOxUZGakBAwbo8ccf14UXXujXeY6dOCyXK7h7ai5u+YPTIVhWcmSb0yEERP6qPk6HEBBx6cH/9agP3xeSdOzAIadDsMxVL3rDjaCZVhksHP9XsXr1ao0bN07r16/XihUrdOLECV1zzTUqKSlxOjQAAEwL1GvE6zrHKxLLly/3+vill15Sq1attGnTJl1++eUORQUAgDWGYa17IlhmbTieSJyusLBQktSiRYtqP19WVqaysjLPx0VFRbUSFwAA/jAMiytbBkki4XjXxs+53W5NnDhRAwcOVPfu3avdJysrSzExMZ4tISGhlqMEAACV6lQiMW7cOH3xxRd64403atxnypQpKiws9Gy5ubm1GCEAAL5xB2ALBnWmayMjI0Pvvvuu1qxZo/POO6/G/cLDwxUeHl6LkQEA4D+3IbktdG6wRLaPDMPQvffeqyVLlmjVqlXq0KGD0yEBAAAfOZ5IjBs3Tq+99prefvttNW3aVHl5eZKkmJgYRUZGOhwdAADm8BrxWjJ37lwVFhbqyiuvVOvWrT3b4sWLnQ4NAADTKt/+aWULBo5XJIxgmd8CAACqcDyRAACgPjL+8z8rxwcDEgkAAGzgtriyZbB0bTg+RgIAAAQvKhIAANjA6qJSLEgFAEADZhgWx0gEyWQEEgkAAGzQUCoSjJEAAACmUZEAAMAGdG0AAADTDFnrngiONIKuDQAAYAEVCQAAbOA2DIuvEQ+OmgSJBAAANmgoS2TTtQEAAEyjIgEAgA0ayjoS9SaRiA47TyGu4L6dQde973QIljVe6XQEgfG37Tc6HUJA3LfyRadDsGzQdU5HEBgxOy90OgTLCko+dTqEAKi97gK3LI6RoGsDAADUd8H9JzwAAHUUszYAAIBpDWXWBokEAAA2YIwEAADAWVCRAADABg2lIkEiAQCADRrKGAm6NgAAgGlUJAAAsIFhsWsjWCoSJBIAANjA7XLL5TK/0LU7SBbJpmsDAACYRkUCAAAbuGXIxawNAABghvGfCaBWjg8GdG0AAADTqEgAAGADt2SxayM4kEgAAGCDhjJrg0QCAAAbuOWWy0IyECyJBGMkAACAaVQkAACwQUOpSJBIAABgA6Z/AgAAnAUVCQAAbMCsDQAAYJoht6VkIFi6NupNIvGLRhersSvM6TAsCRl80ukQLCt4ptzpEAJi4w8up0MIiIIlCU6HYFnL8cH9fV3piqeTnQ7Bsk+iIpwOwTK3cVLfl3zsdBj1Sr1JJAAAqEsMVciwMBTRUEUAo7EPiQQAADY41a1R/8dIMGsDAIB6ZM6cOUpMTFRERIT69eunjRs31rjvjh079Nvf/laJiYlyuVzKzs72+3okEgAA2MDtGW5p9n/+v/Br8eLFyszM1LRp07R582YlJSUpNTVVBQUF1e5/7NgxdezYUTNnzlR8fLyp+ySRAADABqfGSFjb/DVr1iyNHTtW6enp6tatm+bNm6cmTZpowYIF1e5/ySWX6Mknn9SoUaMUHh5u6j4ZIwEAgA0CNUaiqKjIqz08PLzaX/rl5eXatGmTpkyZ4mkLCQlRSkqK1q1bZzqOs6EiAQBAHZaQkKCYmBjPlpWVVe1+hw8fVkVFheLi4rza4+LilJeXZ1t8VCQAALBBoN61kZubq+joaE+72S4Iu5BIAABgA7cqJJlf3M79nzES0dHRXolETWJjYxUaGqr8/Hyv9vz8fNMDKX1B1wYAAPVAWFiY+vTpo5ycHE+b2+1WTk6O+vfvb9t1qUgAAGADJ14jnpmZqbS0NCUnJ6tv377Kzs5WSUmJ0tPTJUljxoxR27ZtPeMsysvL9eWXX3r++8CBA9q6dauioqLUuXNnn65JIgEAgA3chsWuDcP/6Z8jR47UoUOHNHXqVOXl5alXr15avny5ZwDm/v37FRLy386IgwcPqnfv3p6Pn3rqKT311FO64oortGrVKp+uSSIBAEA9kpGRoYyMjGo/d3pykJiYKMPwf+GrnyORAADABk50bTiBRAIAABucSiTMv8EzWBIJZm0AAADTqEgAAGADw3DLbWGwpWEER0WCRAIAABuc6pqwkEgESdcGiQQAADYwTEzfDOTxtYUxEgAAwDQqEgAA2ODUCAm6NgAAgAmnBkvW/8GWdG0AAADTqEgAAGADK4tRBeL42lJvEoku0YbCQ6ytF+60svbJTodg2fOrQ50OISCWFM1zOoSA6Ln6LqdDsOwPTwfHD9Oz6RYT3D+fJOn/9n/sdAgBUHtfh1PvsLCwRLbFd2DUFro2AACAafWmIgEAQF1iddYFszYAAGjATi0oZb57glkbAACg3qMiAQCADaxWFIKlIkEiAQCADRgjAQAATGsoFYk6M0Zizpw5SkxMVEREhPr166eNGzc6HRIAADiLOpFILF68WJmZmZo2bZo2b96spKQkpaamqqCgwOnQAAAwxZDb8hYM6kQiMWvWLI0dO1bp6enq1q2b5s2bpyZNmmjBggVOhwYAgCmGUWF5CwaOJxLl5eXatGmTUlJSPG0hISFKSUnRunXrquxfVlamoqIirw0AADjD8UTi8OHDqqioUFxcnFd7XFyc8vLyquyflZWlmJgYz5aQkFBboQIA4IfKd22Y3XjXhi2mTJmiwsJCz5abm+t0SAAAVGEYbstbMHB8+mdsbKxCQ0OVn5/v1Z6fn6/4+Pgq+4eHhys8PLy2wgMAAGfgeEUiLCxMffr0UU5OjqfN7XYrJydH/fv3dzAyAADMayizNhyvSEhSZmam0tLSlJycrL59+yo7O1slJSVKT093OjQAAExyS3JZOD44xkjUiURi5MiROnTokKZOnaq8vDz16tVLy5cvrzIAEwAA1C11IpGQpIyMDGVkZDgdBgAAgWFYrEgYVCQAAGiwDItdGwZdGwAANGQNY4yE47M2AABA8KIiAQCALQyLRYXgqEiQSAAAYAuroxxIJGqF8Z9RrWXucocjsa6oODje9HYmZe5Qp0MIkOD4Bj4bvi/qjvrxvVEfvi9O3YNRazMi6sMzOzOXUXtP0xbfffcdL+4CAPglNzdX5513ni3nLi0tVYcOHap98aS/4uPjtXfvXkVERAQgMnsEfSLhdrt18OBBNW3aVC6XldGxNSsqKlJCQoJyc3MVHR1tyzVqQ324j/pwD1L9uI/6cA8S91GX1MY9GIah4uJitWnTRiEh9s03KC0tVXm59YpgWFhYnU4ipHrQtRESEmJbVnm66OjooP0G/bn6cB/14R6k+nEf9eEeJO6jLrH7HmJiYmw7d6WIiIg6nwAECtM/AQCAaSQSAADANBIJH4SHh2vatGkKDw93OhRL6sN91Id7kOrHfdSHe5C4j7qkPtxDQxT0gy0BAIBzqEgAAADTSCQAAIBpJBIAAMA0EgkAAGAaicRZzJkzR4mJiYqIiFC/fv20ceNGp0Py25o1azR06FC1adNGLpdL//jHP5wOyW9ZWVm65JJL1LRpU7Vq1UrXX3+9du3a5XRYfpk7d6569uzpWWynf//+WrZsmdNhWTZz5ky5XC5NnDjR6VD88vDDD8vlcnltXbp0cTosvx04cEA33XSTWrZsqcjISPXo0UOfffaZ02H5JTExscrXwuVyady4cU6HBh+QSJzB4sWLlZmZqWnTpmnz5s1KSkpSamqqCgoKnA7NLyUlJUpKStKcOXOcDsW01atXa9y4cVq/fr1WrFihEydO6JprrlFJSYnTofnsvPPO08yZM7Vp0yZ99tlnuvrqq3Xddddpx44dTodm2qeffqrnn39ePXv2dDoUUy666CJ9//33nm3t2rVOh+SXn376SQMHDlTjxo21bNkyffnll3r66afVvHlzp0Pzy6effur1dVixYoUkacSIEQ5HBp8YqFHfvn2NcePGeT6uqKgw2rRpY2RlZTkYlTWSjCVLljgdhmUFBQWGJGP16tVOh2JJ8+bNjb/85S9Oh2FKcXGxcf755xsrVqwwrrjiCmPChAlOh+SXadOmGUlJSU6HYcnkyZONX/ziF06HEXATJkwwOnXqZLjdbqdDgQ+oSNSgvLxcmzZtUkpKiqctJCREKSkpWrdunYORQZIKCwslSS1atHA4EnMqKir0xhtvqKSkRP3793c6HFPGjRuna6+91ut7JNh89dVXatOmjTp27KjRo0dr//79Tofkl6VLlyo5OVkjRoxQq1at1Lt3b73wwgtOh2VJeXm5XnnlFd166622vYgRgUUiUYPDhw+roqJCcXFxXu1xcXEBeTUszHO73Zo4caIGDhyo7t27Ox2OX7Zv366oqCiFh4frrrvu0pIlS9StWzenw/LbG2+8oc2bNysrK8vpUEzr16+fXnrpJS1fvlxz587V3r17ddlll6m4uNjp0Hz2zTffaO7cuTr//PP1/vvv6+6779b48eO1aNEip0Mz7R//+IeOHDmiW265xelQ4KOgf/snGp5x48bpiy++CLr+bEm68MILtXXrVhUWFuqtt95SWlqaVq9eHVTJRG5uriZMmKAVK1YE9dsNBw8e7Pnvnj17ql+/fmrfvr3efPNN3XbbbQ5G5ju3263k5GTNmDFDktS7d2998cUXmjdvntLS0hyOzpwXX3xRgwcPVps2bZwOBT6iIlGD2NhYhYaGKj8/36s9Pz9f8fHxDkWFjIwMvfvuu1q5cmWtvT4+kMLCwtS5c2f16dNHWVlZSkpK0uzZs50Oyy+bNm1SQUGBLr74YjVq1EiNGjXS6tWr9cwzz6hRo0aqqKhwOkRTmjVrpgsuuEB79uxxOhSftW7dukoS2rVr16Droqn07bff6sMPP9Ttt9/udCjwA4lEDcLCwtSnTx/l5OR42txut3JycoK2TzuYGYahjIwMLVmyRB999JE6dOjgdEgB4Xa7VVZW5nQYfhk0aJC2b9+urVu3erbk5GSNHj1aW7duVWhoqNMhmnL06FF9/fXXat26tdOh+GzgwIFVpkHv3r1b7du3dygiaxYuXKhWrVrp2muvdToU+IGujTPIzMxUWlqakpOT1bdvX2VnZ6ukpETp6elOh+aXo0ePev2VtXfvXm3dulUtWrRQu3btHIzMd+PGjdNrr72mt99+W02bNvWMU4mJiVFkZKTD0flmypQpGjx4sNq1a6fi4mK99tprWrVqld5//32nQ/NL06ZNq4xNOeecc9SyZcugGrMyadIkDR06VO3bt9fBgwc1bdo0hYaG6sYbb3Q6NJ/9/ve/14ABAzRjxgzdcMMN2rhxo+bPn6/58+c7HZrf3G63Fi5cqLS0NDVqxK+moOL0tJG67tlnnzXatWtnhIWFGX379jXWr1/vdEh+W7lypSGpypaWluZ0aD6rLn5JxsKFC50OzWe33nqr0b59eyMsLMw499xzjUGDBhkffPCB02EFRDBO/xw5cqTRunVrIywszGjbtq0xcuRIY8+ePU6H5bd33nnH6N69uxEeHm506dLFmD9/vtMhmfL+++8bkoxdu3Y5HQr8xGvEAQCAaYyRAAAAppFIAAAA00gkAACAaSQSAADANBIJAABgGokEAAAwjUQCAACYRiIBAABMI5EAAACmkUgAQerKK6/UxIkTbTn3fffdp6FDh9pybgD1C0tkA0Hqxx9/VOPGjdW0aVNJpxKLXr16KTs72/K5jxw5otDQUM+5AaAmvGINCFItWrSw7dzNmjWz7dwA6he6NoA67K233lKPHj0UGRmpli1bKiUlRSUlJZK8uzZuueUWrV69WrNnz5bL5ZLL5dK+ffvkdruVlZWlDh06KDIyUklJSXrrrbfOeM3Dhw/L5XLpiy++sPv2ANQDVCSAOur777/XjTfeqCeeeEL/8z//o+LiYn388ceqrjdy9uzZ2r17t7p3767p06dLks4991xlZWXplVde0bx583T++edrzZo1uummm3TuuefqiiuuqPa627ZtU3h4uLp06WLr/QGoH0gkgDrq+++/18mTJ/Wb3/xG7du3lyT16NGj2n1jYmIUFhamJk2aKD4+XpJUVlamGTNm6MMPP1T//v0lSR07dtTatWv1/PPPnzGRuOiii9SoET8eAJwdPymAOiopKUmDBg1Sjx49lJqaqmuuuUbDhw9X8+bNfTp+z549OnbsmH75y196tZeXl6t37941Hrd161b16tXLSugAGhASCaCOCg0N1YoVK/TJJ5/ogw8+0LPPPqsHHnhAGzZsUIcOHc56/NGjRyVJ7733ntq2bev1ufDw8BqP27Ztm2677TZrwQNoMBhsCdRhLpdLAwcO1COPPKItW7YoLCxMS5YsqXbfsLAwVVRUeD7u1q2bwsPDtX//fnXu3NlrS0hIqPYc5eXl+ve//62kpCRb7gdA/UNFAqijNmzYoJycHF1zzTVq1aqVNmzYoEOHDqlr167V7p+YmKgNGzZo3759ioqKUosWLTRp0iT9/ve/l9vt1i9+8QsVFhbqX//6l6Kjo5WWllblHP/+97914sQJEgkAPiORAOqo6OhorVmzRtnZ2SoqKlL79u319NNPa/DgwdXuP2nSJKWlpalbt246fvy49u7dq0cffdQze+Obb75Rs2bNdPHFF+v//b//V+05tm7dqvbt27OOBACfsbIlAI+MjAwVFBTozTffdDoUAEGCMRIAVFpaqk2bNulvf/ubUlNTnQ4HQBAhkQCg7OxspaSk6LrrrtOYMWOcDgdAEKFrAwAAmEZFAgAAmEYiAQAATCORAAAAppFIAAAA00gkAACAaSQSAADANBIJAABgGokEAAAwjUQCAACYRiIBAABMI5EAAACm/X9vKkWResHwGgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0, N_s-1, N_s)\n",
    "y = np.linspace(0, (N_t - 1) * dt, N_t)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "plt.pcolormesh(X, Y, occs, cmap='inferno')\n",
    "plt.xlabel(r'site $j$')\n",
    "plt.ylabel(r'time $t$')\n",
    "title = r'(a) ED $\\langle N_j (t) \\rangle$: $N_s = $' + str(N_s)\n",
    "title += r', $J = $' + str(J)\n",
    "# title += r', $h = $' + str(h)\n",
    "title += r', $g = $' + str(g)\n",
    "plt.title(title)\n",
    "plt.colorbar()\n",
    "# plt.savefig('qc-occ-nums.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f233f01d810>"
      ]
     },
     "execution_count": 153,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhsAAAHNCAYAAAC3qLkSAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAUN9JREFUeJzt3XtcFOX+B/DPclkuCqggN0VQM1FB8EpoaSVHNNOoI6BpInnsckAtfprlMbVM0ToaliRipZ5SsfKSlalEaXkgTcDUyltpkApIKSgo6O7z+8PD0sqCuzs7zrJ93r3m9crZZ2a+MyzLd7/PM8+ohBACRERERDKxUzoAIiIism1MNoiIiEhWTDaIiIhIVkw2iIiISFZMNoiIiEhWTDaIiIhIVkw2iIiISFZMNoiIiEhWTDaIiIhIVkw2iIiISFZMNoiIiEhWsiUbr776KoKDg6HVao3eJiMjAx06dEBNTY3F45k3bx5UKhXKy8stvu+mrFmzBiqVCqdPn76tx7U0Q+fx3XffYcCAAWjRogVUKhUOHjzY5HoiIvprkiXZqKysxOLFizFz5kzY2Rl/iIkTJ6K2thYrV640+ZhvvfUWVCoVIiIiTN7WWuXm5mLevHm4ePGiRfdblzjULc7OzvD390d0dDTeeOMNXLp06Zb7uHbtGmJjY/HHH3/g9ddfx3vvvYfAwMBG11NDO3fu1P0Mfvzxxwavjxw5Eu3bt1cgMuDEiRMYM2YM2rdvD1dXVwQHB+Pll19GdXW1pP1+/PHHUKlU+PDDDy0UaUM1NTWYOXMm/P394eLigoiICGRnZ5u1rwULFkClUiEkJMTCUd4+Uq7HDz/8gNjYWHTq1Amurq7w8vLCoEGD8Mknn1gkNmM+t4UQePnll/HNN99Y5JhNscR7p6CgAKNGjUKbNm3g6uqKkJAQvPHGGzJF3IwIGbz++uvC3d1dXLlyxeRtn3vuOREYGCi0Wq1J2w0YMEAEBQUJAOLEiRMNXp87d64AIM6fP29yTFKsXr1aABCnTp0yedvXXnvN7G2Nienll18W7733nnj33XfFwoULxdChQ4VKpRKBgYHi+++/19vm+vXr4sqVK7qfy08//SQAiFWrVum1a2w9NbR48WIBQNjZ2YnFixc3eD0gIEAMHz78tsdVVFQkWrVqJQIDA0VqaqpYuXKlmDhxogAgRo0aJWnf8+fPFwDE0aNHLRRtQ2PGjBEODg5i+vTpYuXKlSIyMlI4ODiIb775xqT9FBcXC1dXV9GiRQvRo0cPmaKVn5Tr8dlnn4no6Ggxb948kZmZKdLS0sQ999wjAIiVK1dKjm3AgAFCrVY3+rkthBBHjx4VAMTGjRslH+9WpL53du7cKdRqtYiIiBBLly4VmZmZYubMmWLGjBkyR279ZEk2evbsKcaPH2/WtgcOHBAARE5OjtHb/PLLLwKA2Lx5s2jbtq2YN29egzZMNhrG9N133zV4LScnR7i4uIjAwEBRXV3d6D727NkjAIgPP/zQqPVSXL582WL7sibjxo0T7u7uYvjw4eLuu+/We+2PP/4QAMTzzz9/2+NasGCBACCOHDmit37ChAkCgPjjjz/M3ndcXJxwcXER169flxqmQfv27RMAxGuvvaZbd+XKFdG5c2cRGRlp0r7i4+PF/fffLwYPHtxskw1LXo86169fF2FhYaJr166SYqv73H755ZeFo6Ojwc9tIYRYt25dk8mIpUi9VhUVFcLHx0c8/PDDQqPRyBlqs2TxZKPuDbRmzZoGr/3222/i8ccfF35+fkKtVougoCDx1FNPiZqaGr12bdq0EVOnTjX6mPPnzxetW7cWNTU14umnnxZdunRp0KYu2fjpp59EbGyscHNz0x3n5gpMZWWlmDZtmggMDBRqtVq0bdtWREVFifz8fF2bgoICMWzYMOHm5iZatGgh7r//fpGXl9fguDcnGwkJCSIwMLDR+G7+981L3X5+++03kZiYKLy9vYVarRbdu3cX77zzjlHXq6lkQwghFi5cKACIzMxMg+eRkJDQIK7Bgwc3ur6OMTHXnfcPP/wgxo4dK1q1aiXCw8PN2seJEydEQkKC8PDwEO7u7mLixImiqqqqwfka87409nr/9NNP4tdff236B/A/ISEhYuDAgSI9PV3Y29uL8vJy3WtfffWVACDWr19v1L4saebMmQYT85kzZwo7OztJyV+3bt1Enz59pIbYqBkzZgh7e3tRUVGht77uPV1UVGTUfvbs2SPs7e3FoUOHLJZsrFmzRoSGhgpnZ2fRu3dvkZubK2JjY0XPnj0l77sxlroeN3vwwQeFj4+PpNjmz58v7O3tRUlJiXjggQcMfm7369evwWeKh4eHpOM2Ruq1WrFihQAgfvzxRyHEjS9JTDrqOVikL+ZPcnNzAQC9e/fWW3/27Fn0798fFy9exBNPPIHg4GCcOXMGH330Eaqrq6FWq3Vte/fujf/+979GH3PdunV45JFHoFarMXbsWKxYsQLfffcd+vXr16BtXFwcgoKCkJqaim+//RZvvPEGLly4gP/85z+6Nk899RQ++ugjJCcno3v37vj999+xd+9e/PTTT+jduzd++OEH3HPPPXB3d8dzzz0HR0dHrFy5Evfeey/27NljkXEjjzzyCI4fP44NGzbg9ddfh5eXFwCgbdu2KC0txV133QWVSoXk5GS0bdsWn3/+OSZNmoTKyko888wzko792GOPYdasWdi1axcmT57c4PUnn3wS7dq1w8KFCzF16lT069cPPj4+aNmypcH1AEyOOTY2Fl26dMHChQshhDBrH3FxcejYsSNSU1NRUFCAt99+G97e3li8eLGujTHvS1OO261bNwwePBi7d+9u8hrX1tbi2LFjmDx5Mh588EEkJSVh+/bteOyxxwAAhw4dAgCEhYUZ8yPTuXbtGioqKoxq26ZNG4Njqu69914sXrwYkyZNwksvvQRPT0/k5uZixYoVmDp1Klq0aGFSTHVqampw4sQJjB8/XrbYCwsLceedd8Ld3V1vff/+/QEABw8eREBAQJP71mg0mDJlCv7xj38gNDTUqHhuJTU1FbNmzcLDDz+MKVOm4ODBgxg1ahQ8PDxw1113GdzGWq4HAFRVVeHKlSuoqKjAtm3b8PnnnyM+Pt6o2Bqzbt06DBo0CD4+PoiLi8PEiRMbfG7PnDkT8+bNQ01NDebMmQMAaNWqVYN9WcO1+uKLL+Du7o4zZ84gJiYGx48fR4sWLfDYY4/h9ddfh7Ozs1Hx2SxLZy+zZ88WAMSlS5f01k+YMEHY2dkZ/DZ98/iMJ554Qri4uBh1vLpul+zsbN2+2rdvL6ZNm6bXru7b7s19zv/85z8FAL0xCh4eHiIpKanRY8bExAi1Wi1+/vln3bqzZ88KNzc3MWjQIL225lY2hGi8G2XSpEnCz89P75uwEDf6Gz08PJrs/vhzTI1VNoS4cQ169erV6HnUffO+ubuksfXGxlx3HcaOHdsgJlP38fjjj+u1e/jhh4Wnp6feOmPel6Zcb9xUzWlMYWGhACAyMjKEEEKEhoaK2NhY3euPP/64cHJyMrm7oe76G7M01T03f/584eLiotf+X//6l0mx3KzunJcuXSpb7D169BD3339/g/U//PCD3vVuyvLly4WHh4coKysTQgjJlY38/Hxhb28vXnjhBb31kyZNEgBEamqqwe2s5XoIIcSTTz6pO5adnZ0YPXq0pO60us/tuuNfvHhRqNXqBp/bQgjRoUMHMXHixCb3Zw3XqmfPnsLV1VW4urqKKVOmiE2bNokpU6YIAGLMmDFNbvtXYPHKxu+//w4HBwe0bNlSt06r1WLr1q0YOXIk+vbt22AblUql9+/WrVvjypUrqK6uhqura5PHW7duHXx8fHDffffp9hUfH4/3338fS5Ysgb29vV77pKQkvX9PmTIFb731FrZv346ePXsCuJE579u3D2fPnoW/v79ee41Gg127diEmJgadOnXSrffz88Ojjz6KVatWobKyskF2bClCCGzatAlxcXEQQujdyhsdHY2srCwUFBRg4MCBko7TsmVLo+5KMYY5MT/11FMW38c999yDLVu26H4+xrwvTT2u+F8V5lbqKhd177kHH3wQ6enpuHbtGhwdHfH999+jR48eDd6/txIWFmb06HlfX99GXwsKCsKgQYPw97//HZ6envjss8+wcOFC+Pr6Ijk52aSY6hw+fBgAGq0WWCL2K1euwMnJqcH6um+VV65caXK/v//+O+bMmYMXX3wRbdu2NSqWW1m4cCH8/Pzw0ksv6a2/55578M4771j19ajzzDPPYPTo0Th79iw++OADaDQa1NbWGrWtIevWrYODgwP+/ve/AwA8PDwwbNgwZGVl6X1uV1RUoKioSPd70hhruFaXL19GdXU1nnrqKd3dJ4888ojuDsuXX34ZXbp0MSpGW2TxZMOQ8+fPo7Ky0ujbx+o+sG9OQm6m0WiQlZWF++67D6dOndKtj4iIwJIlS5CTk4OhQ4fqbXPzD7tz586ws7PTmz/i1VdfRUJCAgICAtCnTx888MADmDBhAjp16oTz58+juroaXbt2bRBPt27doNVqUVxcjB49ehh1rqY6f/48Ll68iMzMTGRmZhpsU1ZWJvk4ly9fhre3t+T9AObF3LFjR8n76NChg96/W7duDQC4cOEC3N3djXpfynW9v//+e6hUKt0fmgcffBCpqan4+uuvce+99+KHH37AmDFjTN5v69atERUVZfJ2f5aVlYUnnngCx48f1916+8gjj0Cr1WLmzJkYO3YsPD09Td5vXbLR2B8OS8Tu4uJicJ6eq1ev6l5vyuzZs9GmTRtMmTJFUhx1ampq8Pnnn+Ppp5+Go6Oj3muXL18G0HjyZQ3Xo05wcDCCg4MBABMmTMDQoUMxcuRI7Nu375af0zer+9y+//77dd3DABAfH49t27bpfW7fnJQ3xhquVd3rY8eO1Vv/6KOPYuXKlcjLy2OyYUmenp64fv06Ll26BDc3N7P2ceHCBbi6ut7yh/vll1/i3LlzyMrKQlZWVoPX161b1yDZuJmhX5S4uDjdt+Bdu3bhtddew+LFi7F582b06tXLtJMx4njAjV9AY9RNkjZ+/HgkJCQYbHOrX8xb+e2331BRUYE77rhD0n7qmBPzzT97c/bRWFXA2OqDucc1xqFDh9CpUyddBfCuu+6Cl5cXPvnkE/j7++Pq1asmj9cAbowF+eOPP4xq27ZtW4PX6K233kKvXr0azPExatQorFmzBoWFhWZ9sB8+fBje3t6NJrGWiN3Pzw9nzpxpsP7cuXMA0KBS+WcnTpxAZmYm0tLScPbsWd36q1ev4tq1azh9+jTc3d3Rpk0bo2IEgJ9//hnV1dUGf5a//PIL3N3dGyTFdZS+Hk0ZPXo0nnzySRw/ftzgF6+m1H1uv/LKK3rrR40aBRcXF73PbWPHLlnDtfL398cPP/ygG6dWp+79fuHCBaPis1UWTzbqst9Tp07pPoTbtm0Ld3d3HDlyxKh9nDp1Ct26dbtlu3Xr1sHb2xvp6ekNXtu8eTO2bNmCjIwMvT9cJ06c0PvWfPLkSWi1WgQFBelt7+fnh3/+85/45z//ibKyMvTu3RsLFizAnj174OrqimPHjjU45tGjR2FnZ9fkIKLWrVsbnKTr119/bbDOUGLStm1buLm5QaPRSM7kG/Pee+8BuNFNYAmWiFmO8zbmfSnX9T506JBel4+dnR2GDx+OTz75RDdg8M9JzNGjR5GcnIyCggIIIfDYY48ZnCgoNzdX16V4K6dOnWrwvgduDMStqwL92bVr1wAA169fN2r/Nzt8+HCTiZklYg8PD8dXX33VoCtz3759utcbc+bMGWi1WkydOhVTp05t8HrHjh0xbdo0pKWlGRUjUF96v3lAYk1NDTZs2NBkVU3p69GUuvMydlDmn61btw6Ojo54+OGH9da3bNkSDzzwgN7n9qFDh+Dn56dXATHEGq5Vnz59kJ2djTNnzuglYHWJq6W65ZoriycbkZGRAIADBw7oPljs7OwQExOD999/HwcOHGjQPy6E0PvDWlBQgHHjxjV5nCtXrmDz5s2IjY3F6NGjG7zu7++PDRs2YNu2bXqjptPT0/WqHW+++SYAYPjw4QBuVBguX74MDw8PXRtvb2/4+/ujpqYG9vb2GDp0KD7++GOcPn1a96YtLS3F+vXrcffddzc5XqNz586oqKjAoUOHdNfn3Llz2LJlS4O2daP+/5yc2Nvb4+9//zvWr1+PI0eONPiwOn/+vKQ39Zdffon58+ejY8eOt/wZGMsSMctx3sa8L0097tGjR+Hq6trot1UAKCkpQVlZWYNvaw8++CDee+89bNiwAYD+t7lx48Zh5syZiI2NxaVLl3DixAmD+7ZE3/Wdd96JXbt24fjx47jzzjt16zds2AA7OzuzKjkXLlzAmTNnEBcX12gbS8Q+evRo/Pvf/0ZmZiamT58O4MYf9tWrVyMiIkL3RaC6uhpFRUXw8vLS/SELCQkx+Hs4e/ZsXLp0CcuWLUPnzp2Niq9O3ey5e/fu1ft9WrRoEc6dO4eRI0c2uq3S1wO40UV4cyXq2rVr+M9//gMXFxd0797dqPjq1H1u/+1vfzOY0MbFxWHTpk26z+2ioiKjZtG9ndcKMHy94uLisGjRIrzzzju4//77dW3ffvttODg44N577zUqPpslx6jTkJCQBncT/Pbbb8LX11e4urqKZ555RqxcuVLMmzdP9OjRQ1y4cEHXrm6U8hdffNHkMbKysgQAsXXrVoOvazQa0bZtWzFy5EghRP0dCqGhoWLkyJEiPT1djB8/XgAQjz76qG67CxcuiBYtWoiEhATdDHBxcXECgFiyZIkQQogjR46IFi1aiHbt2okFCxaIxYsXi06dOgknJyfx7bff6sVx810c5eXlokWLFqJTp04iLS1NLFy4UAQEBIjevXs3uBtl//79AoB44IEHxH/+8x+xYcMGcfnyZVFSUiICAwOFq6urmDZtmli5cqVITU0VsbGxonXr1k1etz/HVDeD6OrVq8WiRYt0M4gGBQWJw4cPN3kept6NYmzMTU2+JnUfhiZYM+Z9acr1hhF3o+zYsUMANyah+7OLFy8KR0dHoVKphL+/v95rHh4e4r333pNtMqw/q5tjwtvbW7z88ssiPT1dDB8+XAAQ//jHPxq0N+ac6yZ7W716tTxB/0lsbKxwcHAQM2bMECtXrhQDBgwQDg4OYs+ePbo2de/TuXPn3nJ/jd2NYsx5CyHEsGHDhJ2dnUhOTtZ9nnTu3FkAEG+++aYpp2YWKdcjJiZG3H///WLevHli1apVYv78+SI4OFjv87COMdej7nP7gQceEKmpqQ2WF198UQDQfW4//fTTQq1Wi8WLF4v33ntPHDhwwCLXpDHGXCshGr9ejz/+uAAg4uLiRHp6uoiNjRUAGtyJ9FckS7KxdOlS0bJlywa3YP76669iwoQJom3btsLJyUl06tRJJCUl6U2eNHPmTNGhQ4dbTlc+cuRI4ezsbHCSpjoTJ04Ujo6Oory8XPcH6McffxSjR48Wbm5uonXr1iI5OVlvUq+amhoxY8YMERYWppuwKywsTLz11lt6+y4oKBDR0dGiZcuWwtXVVdx3330iNze3QQyG/sDt2rVLhISECLVaLbp27Sref/99g7e+CnHjFsR27doJOzs7vf2UlpaKpKQkERAQIBwdHYWvr68YMmSI3kRcjamLqW5Rq9XC19dX/O1vfxPLli0TlZWVtzwPU5MNY2O+1UyvUvbR2Gyuxrwvjb3exnzgvvrqqwKAOHnyZIPX7rvvPgFADBs2TG/99u3bxcCBA4WPj4+YMWOGuHbtWpPHkGrfvn1i+PDhwtfXVzg6Ooo777xTLFiwoMFxL126ZNStfcuXLxcAREFBgZxhCyFuzPo4ffp04evrK5ycnES/fv3Ejh079NpITTaMPW8hbrx3YmJihLu7u/Dz8xPTpk0TW7ZsEQAafDmRg5TrsWHDBhEVFSV8fHyEg4ODaN26tYiKihIff/yxXjtjr8fIkSONuj217nP7zJkzus9ZAOKNN96wyDVpjDHXSojGr1dtba2YN2+eCAwMFI6OjuKOO+4Qr7/+uqwxNxeyJBsXL14Ubdq0EW+//bZJ2129elX4+vqKtLQ0OcIiavZOnz4tOnToIHbu3Kl0KEKIG8/OUKlU4tChQ022e/TRR4Wbm5vsSdLtYux5NyYpKUl4eXnxetBfhixPffXw8MBzzz2H1157zaRHzK9evRqOjo4N5kcg+ivbtGmT7tbuCxcuoLa2VjcQW2lfffUVxowZ0+RMm+fPn8e2bdswevRoODjclrvtZWfMeQM3+vzFTXc/ff3111i5ciWeeuqpv9z1oL8ulbj5N4GIrMrUqVPxwQcf4PLly+jcuTMWLlyIESNGKB3WLR0+fBg5OTlYsWIFSkpKcPDgwQbzp9i6bdu24eWXX8bo0aPRqlUrFBQUYPXq1ejVqxf27Nlj9DwXRM0dkw0ikkVCQgI2b96MgQMHYvHixWbNG9Lc7d27F8899xx+/PFHVFdXIygoCHFxcXjhhRfMfsYMUXPEZIOIiIhkJcuYDSIiIqI6TDaIiIhIVs1+KLRWq8XZs2fh5uZm8gOBiIjor0UIgUuXLsHf37/BNPKWdPXqVUlPxq2jVqt1T55tzpp9snH27Nkmn0VCRER0s+LiYqOmQjfH1atX0bGjL0pKTH92zM18fX1x6tSpZp9wNPtko+7Jsv8JGwlXe8dbtLZuP19spXQIkj0YXqh0CBZx9+fGzw9jzS5ePal0CJL1cLb+23yNseMfe5UOQbKMj5v/z6JGW4vXfltr9lPJjVFbW4uSkgqcLl4Gd3fzb2+urLyCoIBpqK2tZbKhtLquE1d7R7Ro5smGi71a6RAkc3M0/Fj35sZWuuRs4TzsVc3797qOu1PzHyLnbNf8P6Pq3I7fjZYtndCypZPZ25syKaa1a/bJBhERkTUS4jqEuC5pe1vBZIOIiEgGQmgghEbS9rai+df1iIiIyKqxskFERCQDrbgOrYSuECnbWhsmG0RERDLgmI167EYhIiIiWbGyQUREJIMbA0SlVDZsZ4Aokw0iIiIZCO11CK2EZEPCttaG3ShEREQkK1Y2iIiI5CCu31ikbG8jmGwQERHJgHej1GM3ChEREcmKlQ0iIiI5aK8D2mvStrcRTDaIiIhkcKMbxfwnYdtSNwqTDSIiIjlorwNa85MNW6pscMwGERERyYqVDSIiIjmwsqHDZIOIiEgWGolzZdjOdOXsRiEiIiJZsbJBREQkA5X2OlRa87/Tq9iNQkRERE3SXgckJBu2NGaD3ShEREQkK8UrGxqNBvPmzcP777+PkpIS+Pv7Y+LEiZg9ezZUKpXR+/lvSVs42alljFR+ub/XKB2CZD9UDFA6BIsY5650BJax0y5Q6RAk2/bAcaVDsAg7x+b/LVVtJ5QOQTItbuM5sLKho3iysXjxYqxYsQJr165Fjx49cODAASQmJsLDwwNTp05VOjwiIiKzqMR1qISEMRucQdRycnNz8dBDD2HEiBEAgKCgIGzYsAH79+9XODIiIiKyBMXHbAwYMAA5OTk4fvxGqfT777/H3r17MXz4cIPta2pqUFlZqbcQERFZHa0W0GokLFqlz8BiFK9sPP/886isrERwcDDs7e2h0WiwYMECjBs3zmD71NRUvPTSS7c5SiIiItPcuPXV+LGHhra3FYpXNj744AOsW7cO69evR0FBAdauXYt///vfWLt2rcH2L7zwAioqKnRLcXHxbY6YiIjICJKqGv9bbITilY0ZM2bg+eefx5gxYwAAoaGh+PXXX5GamoqEhIQG7Z2cnODk5HS7wyQiIiIzKZ5sVFdXw85Ov8Bib28PrQ31VRER0V+Q9jogoRuFt75a0MiRI7FgwQJ06NABPXr0QGFhIZYuXYrHH39c6dCIiIjMptJqJE5Xzm4Ui3nzzTfx4osv4p///CfKysrg7++PJ598EnPmzFE6NCIiIrIAxZMNNzc3pKWlIS0tTelQiIiILEdopM0gKljZICIioiaotFpJXSEqGxq7qPitr0RERGTbWNkgIiKSg1Yj8W4UdqMQERFRE27cjSJlBlHbSTbYjUJERESyYrJBREQkB4WmK09PT0dQUBCcnZ0RERFxy6eof/jhhwgODoazszNCQ0Oxffv2Bm1++uknjBo1Ch4eHmjRogX69euHoqIio2NiskFERCSDG90o0hZTbdy4ESkpKZg7dy4KCgoQFhaG6OholJWVGWyfm5uLsWPHYtKkSSgsLERMTAxiYmJw5MgRXZuff/4Zd999N4KDg7F7924cOnQIL774IpydnY2/FkIIYfLZWJHKykp4eHjgWf8n4WSnVjocSXJ/r1E6BMm6uRn/5rNmahtJw3debv4PKsx+4DelQ7AIz87N/2fx1sZHlA5BsqvaWsw5/TYqKirg7u4uyzHq/i6V7ekD95bmD42svHwd3oPzTYo1IiIC/fr1w/LlywEAWq0WAQEBmDJlCp5//vkG7ePj41FVVYVPP/1Ut+6uu+5CeHg4MjIyAABjxoyBo6Mj3nvvPbPPxUY+UomIiGxTZWWl3lJTY/iLaW1tLfLz8xEVFaVbZ2dnh6ioKOTl5RncJi8vT689AERHR+vaa7VafPbZZ7jzzjsRHR0Nb29vREREYOvWrSadA5MNIiIiGai04n8Te5m73Oh4CAgIgIeHh25JTU01eLzy8nJoNBr4+Pjorffx8UFJSYnBbUpKSppsX1ZWhsuXL2PRokUYNmwYdu3ahYcffhiPPPII9uzZY/S14K2vREREctBqACmTgP5vzEZxcbFeN4qTk5PEwEwI4X+zmD700EN49tlnAQDh4eHIzc1FRkYGBg8ebNR+mGwQERFZMXd3d6PGbHh5ecHe3h6lpaV660tLS+Hr62twG19f3ybbe3l5wcHBAd27d9dr061bN+zdu9foc2A3ChERkRyExNteTXwQm1qtRp8+fZCTk6Nbp9VqkZOTg8jISIPbREZG6rUHgOzsbF17tVqNfv364dixY3ptjh8/jsDAQKNjs5nKRkb5TqhUzTt3Kvhbe6VDkCz5y+63btQMhHo07zub6pyszlY6BMle+epxpUOwjK8Mf9g3J/urflc6BMk04tptO5ZKaKESEmYQFab3waSkpCAhIQF9+/ZF//79kZaWhqqqKiQmJgIAJkyYgHbt2unGfUybNg2DBw/GkiVLMGLECGRlZeHAgQPIzMzU7XPGjBmIj4/HoEGDcN9992HHjh345JNPsHv3bqPjsplkg4iI6K8uPj4e58+fx5w5c1BSUoLw8HDs2LFDNwi0qKgIdnb1X8wHDBiA9evXY/bs2Zg1axa6dOmCrVu3IiQkRNfm4YcfRkZGBlJTUzF16lR07doVmzZtwt133210XDYzz4aLOoiVDSvAyoZ1WVb6rtIhSDbJ00YqGzbAViobR658cFvm2Sjf2RXuLezN30+VBl7Rx2SN9XZhZYOIiEgOWq3Ep75KuZXFujTvUgARERFZPVY2iIiI5MDKhg6TDSIiIhncmAVU2va2gskGERGRHLRaiTOI2k6ywTEbREREJCtWNoiIiOTAyoYOkw0iIiI5MNnQYTcKERERyYqVDSIiIjkIDaCVMEm3Gc9GsVZMNoiIiGTAW1/rsRuFiIiIZMXKBhERkRw4QFSHyQYREZEcmGzosBuFiIiIZMXKBhERkRy0Qlp1QsqdLFaGyQYREZEctEJiNwqTDSIiImqK5EfM206ywTEbREREJCtWNoiIiOTAyoYOkw0iIiI5cMyGDrtRiIiISFY2U9lwc/SFnap5n45aXat0CJL1beOodAgW8Upxf6VDsIjTHs1/UiBfF9v4drfqj0KlQ5Csl6qn0iFIdg21OHK7Dia0gJDQjSJs470P2FCyQUREZFWExG4UG0o22I1CREREsmJlg4iISA4cIKrDZIOIiEgOTDZ02I1CREREsmJlg4iISAZCe2ORsr2tYLJBREQkB3aj6DDZICIikoMWEpMNSwWiPI7ZICIiIlmxskFERCQHVjZ0mGwQERHJQfxvkbK9jWA3ChEREcmKlQ0iIiIZCK0KQmv+g9h46ysRERE1jWM2dNiNQkRERLJiZYOIiEgOQgVI6EaxpQGiTDaIiIhkwDEb9diNQkRERLJiZYOIiEgOWondKDZU2WCyQUREJAehurGYvb3lQlEau1GIiIhkUDdmQ8pijvT0dAQFBcHZ2RkRERHYv39/k+0//PBDBAcHw9nZGaGhodi+fbve6xMnToRKpdJbhg0bZlJMTDaIiIhsxMaNG5GSkoK5c+eioKAAYWFhiI6ORllZmcH2ubm5GDt2LCZNmoTCwkLExMQgJiYGR44c0Ws3bNgwnDt3Trds2LDBpLiYbBAREclBayd9MdHSpUsxefJkJCYmonv37sjIyICrqyveffddg+2XLVuGYcOGYcaMGejWrRvmz5+P3r17Y/ny5XrtnJyc4Ovrq1tat25tUlxMNoiIiORQN0BUymKC2tpa5OfnIyoqSrfOzs4OUVFRyMvLM7hNXl6eXnsAiI6ObtB+9+7d8Pb2RteuXfH000/j999/Nyk2mxkgmuLTHc72aqXDkGT6F45KhyDZwLbXlQ7BIpwdpiodgkXk3nOX0iFIdq6yldIhWMSK368qHYJkO6vXKB2CZEI0v1GXlZWVev92cnKCk5NTg3bl5eXQaDTw8fHRW+/j44OjR48a3HdJSYnB9iUlJbp/Dxs2DI888gg6duyIn3/+GbNmzcLw4cORl5cHe3t7o87BZpINIiIiayKECkLC3Sh1eVFAQIDe+rlz52LevHkSIjPNmDFjdP8fGhqKnj17onPnzti9ezeGDBli1D6YbBAREclBaydxno0b2UZxcTHc3d11qw1VNQDAy8sL9vb2KC0t1VtfWloKX19fg9v4+vqa1B4AOnXqBC8vL5w8edLoZINjNoiIiKyYu7u73tJYsqFWq9GnTx/k5OTo1mm1WuTk5CAyMtLgNpGRkXrtASA7O7vR9gDw22+/4ffff4efn5/R58DKBhERkQyEFhKfjWL6+JKUlBQkJCSgb9++6N+/P9LS0lBVVYXExEQAwIQJE9CuXTukpqYCAKZNm4bBgwdjyZIlGDFiBLKysnDgwAFkZmYCAC5fvoyXXnoJf//73+Hr64uff/4Zzz33HO644w5ER0cbHReTDSIiIjlIfuqr6dvGx8fj/PnzmDNnDkpKShAeHo4dO3boBoEWFRXBzq6+U2PAgAFYv349Zs+ejVmzZqFLly7YunUrQkJCAAD29vY4dOgQ1q5di4sXL8Lf3x9Dhw7F/PnzG62wGMJkg4iIyIYkJycjOTnZ4Gu7d+9usC42NhaxsbEG27u4uGDnzp2SY2KyQUREJAPpd6NIqIpYGasYIHrmzBmMHz8enp6ecHFxQWhoKA4cOKB0WEREROZTYAZRa6V4ZePChQsYOHAg7rvvPnz++edo27YtTpw4YfJUqERERNZEysPU6ra3FYonG4sXL0ZAQABWr16tW9exY0cFIyIiIiJLUrxGs23bNvTt2xexsbHw9vZGr169sGrVqkbb19TUoLKyUm8hIiKyNnVjNqQstkLxZOOXX37BihUr0KVLF+zcuRNPP/00pk6dirVr1xpsn5qaCg8PD91y8zSuREREVoFjNnQUPxOtVovevXtj4cKF6NWrF5544glMnjwZGRkZBtu/8MILqKio0C3FxcW3OWIiIiIyheJjNvz8/NC9e3e9dd26dcOmTZsMtm/saXdERETWhANE6ymebAwcOBDHjh3TW3f8+HEEBgYqFBEREZF0nGejnuLdKM8++yy+/fZbLFy4ECdPnsT69euRmZmJpKQkpUMjIiIiC1A82ejXrx+2bNmCDRs2ICQkBPPnz0daWhrGjRundGhERETm4wBRHcW7UQDgwQcfxIMPPqh0GERERBbDMRv1bCdtIiIiIqtkFZUNIiIiW8MBovWYbBAREclBSBx3ISwXitKYbBAREcmAYzbqccwGERERycpmKhsCzf+hNZ3dmnf8ADD957eVDsEiVnRNVDoEixj038+VDkGyKd4PKx2CRXihRukQJAtwCVY6BMk0ohaHrnxwW44lhLRxF4LdKERERNQkid0oYDcKERERkXFY2SAiIpKBEHYQwvzv9MKG+lGYbBAREclBq5LWFcJuFCIiIiLjsLJBREQkA84gWo/JBhERkQw4qVc9dqMQERGRrFjZICIikgHvRqnHZIOIiEgG7Eapx2SDiIhIBhwgWo9jNoiIiEhWrGwQERHJgJWNekw2iIiIZCCExDEbNpRssBuFiIiIZMXKBhERkQx462s9JhtEREQy4K2v9diNQkRERLJiZYOIiEgGvBulHpMNIiIiGTDZqMdkg4iISAZCK23chdBaMBiFccwGERERyYqVDSIiIhmwG6Uekw0iIiIZSJ9nw3Y6H2wm2Thb7QAnO0elw5BEYwPzt0S4PKZ0CBZRVO2kdAgWMdb9YaVDkOyrCxVKh2ARkS7tlA5Bsk3VOUqHIJkQGqVD+EuynbSJiIjIimiFSvJijvT0dAQFBcHZ2RkRERHYv39/k+0//PBDBAcHw9nZGaGhodi+fXujbZ966imoVCqkpaWZFBOTDSIiIjn8bwZRcxeYcSfLxo0bkZKSgrlz56KgoABhYWGIjo5GWVmZwfa5ubkYO3YsJk2ahMLCQsTExCAmJgZHjhxp0HbLli349ttv4e/vb3JcTDaIiIhsxNKlSzF58mQkJiaie/fuyMjIgKurK959912D7ZctW4Zhw4ZhxowZ6NatG+bPn4/evXtj+fLleu3OnDmDKVOmYN26dXB0NH3IApMNIiIiGdTdjSJlMUVtbS3y8/MRFRWlW2dnZ4eoqCjk5eUZ3CYvL0+vPQBER0frtddqtXjssccwY8YM9OjRw6SY6tjMAFEiIiJrYqlbXysrK/XWOzk5wcmp4SD28vJyaDQa+Pj46K338fHB0aNHDR6jpKTEYPuSkhLdvxcvXgwHBwdMnTrVrPMAWNkgIiKyagEBAfDw8NAtqampt+3Y+fn5WLZsGdasWQOVyvzEiZUNIiIiGViqslFcXAx3d3fdekNVDQDw8vKCvb09SktL9daXlpbC19fX4Da+vr5Ntv/mm29QVlaGDh066F7XaDT4v//7P6SlpeH06dNGnQsrG0RERDLQCjvJCwC4u7vrLY0lG2q1Gn369EFOTv18KFqtFjk5OYiMjDS4TWRkpF57AMjOzta1f+yxx3Do0CEcPHhQt/j7+2PGjBnYuXOn0deClQ0iIiIZCKGS9iA2M6oiKSkpSEhIQN++fdG/f3+kpaWhqqoKiYmJAIAJEyagXbt2uq6YadOmYfDgwViyZAlGjBiBrKwsHDhwAJmZmQAAT09PeHp66h3D0dERvr6+6Nq1q9FxMdkgIiKyEfHx8Th//jzmzJmDkpIShIeHY8eOHbpBoEVFRbCzq+/UGDBgANavX4/Zs2dj1qxZ6NKlC7Zu3YqQkBCLxsVkg4iISAZKPYgtOTkZycnJBl/bvXt3g3WxsbGIjY01ev/GjtP4MyYbREREMuBTX+txgCgRERHJipUNIiIiGUh5mFrd9raCyQYREZEM2I1Sj90oREREJCtWNoiIiGTAykY9syobAwYMaPBgGCIiIqpXN2ZDymIrzEo2vv32W1y9erXB+srKSsycOVNyUERERGQ7TEo2Ro8ejUWLFkGlUqGsrKzB61VVVfj3v/9tseCIiIiaKyHqu1LMW5Q+A8sxacxGhw4d8Omnn0IIgbCwMHh6eiIsLAxhYWEIDw/HsWPH4OfnJ1esREREzQbHbNQzKdlYunQpgBtPlvvvf/+Ls2fPorCwEAcPHsSWLVug1Wrx6quvyhIoERFRcyIkjrv4yyYbdaqqquDo6AgAeOihhywaEBEREdkWs5KNukSDiIiIDGM3Sj2bmWcju/I87FXNOwn66crnSocg2SDnsUqHYBE/X7KNX/LS2hqlQ5Dsot1FpUOwiPv9mv8cikd/vkvpECS7jlp8h6O35VhMNuo1/3c/ERERWTWbqWwQERFZEz6IrR6TDSIiIhmwG6We2d0o33zzDcaPH4/IyEicOXMGAPDee+9h7969FguOiIiImj+zko1NmzYhOjoaLi4uKCwsRE3NjUFoFRUVWLhwoUUDJCIiao74bJR6ZiUbr7zyCjIyMrBq1Sq922AHDhyIgoICiwVHRETUXAmoJC+2wqxk49ixYxg0aFCD9R4eHrh48aLUmIiIiMiGmJVs+Pr64uTJkw3W7927F506dZIcFBERUXMn7SFs0gaXWhuzko3Jkydj2rRp2LdvH1QqFc6ePYt169Zh+vTpePrppy0dIxERUbPDMRv1zLr19fnnn4dWq8WQIUNQXV2NQYMGwcnJCdOnT8eUKVMsHSMREVGzw1tf65mVbKhUKvzrX//CjBkzcPLkSVy+fBndu3dHy5YtLR0fERERNXNmT+p19epVHDp0CGVlZdBqtSgpKdG9NmrUKIsER0RE1FxpIXEGURu6G8WsZGPHjh147LHH8Pvvvzd4TaVSQaPRSA6MiIioOWM3Sj2zBohOmTIFcXFxOHfuHLRard7CRIOIiIj+zKzKRmlpKVJSUuDj42PpeIiIiGyCFipJXSG21I1iVmVj9OjR2L17t4VDARYtWgSVSoVnnnnG4vsmIiK6raTOsWFD3ShmVTaWL1+O2NhYfPPNNwgNDdWbshwApk6davI+v/vuO6xcuRI9e/Y0JyQiIiKyUmYlGxs2bMCuXbvg7OyM3bt3Q6Wqz75UKpXJycbly5cxbtw4rFq1Cq+88oo5IREREVkVqRNz2dKkXmZ1o/zrX//CSy+9hIqKCpw+fRqnTp3SLb/88ovJ+0tKSsKIESMQFRV1y7Y1NTWorKzUW4iIiKwNpyuvZ1Zlo7a2FvHx8bCzMytX0ZOVlYWCggJ89913RrVPTU3FSy+9JPm4REREdHuYlS0kJCRg48aNkg9eXFyMadOmYd26dXB2djZqmxdeeAEVFRW6pbi4WHIcRERElqa1wGIrzKpsaDQavPrqq9i5cyd69uzZYIDo0qVLjdpPfn4+ysrK0Lt3b719f/3111i+fDlqampgb2+vt42TkxOcnJzMCZuIiOi24aRe9cxKNg4fPoxevXoBAI4cOaL32p8Hi97KkCFDcPjwYb11iYmJCA4OxsyZMxskGkRERM2FVkgb5KkVFgxGYWYlG1999ZVFDu7m5oaQkBC9dS1atICnp2eD9bdyqnYfVCrpY0iU9H6PR5QOQbIvztlG1Sm+Y8mtGzUDo29K5pslG6kl55b9TekQJPuu9n2lQ5BMCBv6C96MmP0gNiIiImqcgApCwiygUra1NkYnGykpKZg/fz5atGiBlJSUJtsaO2bDEDlmJiUiIrrdOM9GPaOTjcLCQly7dk33/40xZcwGERER2T6jk40/j9NYu3Yt2rdv32CeDSEEb0UlIiJC3QBRadvbCrNGVHbs2BHl5eUN1v/xxx/o2LGj5KCIiIiau7oxG1IWW2FWstHYaN7Lly8bPTkXERER/TWYdDdK3cBQlUqFOXPmwNXVVfeaRqPBvn37EB4ebtEAiYiImiOlBoimp6fjtddeQ0lJCcLCwvDmm2+if//+jbb/8MMP8eKLL+L06dPo0qULFi9ejAceeED3+rx585CVlYXi4mKo1Wr06dMHCxYsQEREhNExmVTZKCwsRGFhIYQQOHz4sO7fhYWFOHr0KMLCwrBmzRpTdklERGSThJC+mGrjxo1ISUnB3LlzUVBQgLCwMERHR6OsrMxg+9zcXIwdOxaTJk1CYWEhYmJiEBMTozdh55133only5fj8OHD2Lt3L4KCgjB06FCcP3/e6LhUwowZThITE7Fs2TK4u7ubuqnFVVZWwsPDA04O7Zr9pF7vBg9ROgTJvjin/HvCEjipF1naePfmP6nXqt9tY1IvjfYCKioqZPsbVvd3aVW3R+FqrzZ7P9WaWkz+ab1JsUZERKBfv35Yvnw5AECr1SIgIABTpkzB888/36B9fHw8qqqq8Omnn+rW3XXXXQgPD0dGRobBY9Sd3xdffIEhQ4z7u2XWX+fVq1dbRaJBRERkrQRU0EpY6gaIVlZW6i01NTUGj1dbW4v8/HxERUXp1tnZ2SEqKgp5eXkGt8nLy9NrDwDR0dGNtq+trUVmZiY8PDwQFhZm9LVo3qUAIiIiK1X3IDYpCwAEBATAw8NDt6Smpho8Xnl5OTQaDXx8fPTW+/j4oKTEcLW2pKTEqPaffvopWrZsCWdnZ7z++uvIzs6Gl5eX0deC05UTERHJwFIDRIuLi/V6E5R48vl9992HgwcPory8HKtWrUJcXBz27dsHb29vo7ZnZYOIiMiKubu76y2NJRteXl6wt7dHaWmp3vrS0lL4+voa3MbX19eo9i1atMAdd9yBu+66C++88w4cHBzwzjvvGH0OTDaIiIhkICywmKLuttScnBzdOq1Wi5ycHERGRhrcJjIyUq89AGRnZzfa/s/7bWzsiCHsRiEiIpKBEvNspKSkICEhAX379kX//v2RlpaGqqoqJCYmAgAmTJiAdu3a6cZ9TJs2DYMHD8aSJUswYsQIZGVl4cCBA8jMzAQAVFVVYcGCBRg1ahT8/PxQXl6O9PR0nDlzBrGxsUbHxWSDiIjIRsTHx+P8+fOYM2cOSkpKEB4ejh07dugGgRYVFek912zAgAFYv349Zs+ejVmzZqFLly7YunUrQkJCAAD29vY4evQo1q5di/Lycnh6eqJfv3745ptv0KNHD6PjYrJBREQkA+3/FinbmyM5ORnJyckGX9u9e3eDdbGxsY1WKZydnbF582YzI6nHZIOIiEgGf7591dztbQUHiBIREZGsWNkgIiKSgVIPYrNGTDaIiIhkYM7tqzdvbyvYjUJERESyYmWDiIhIBuxGqcdkg4iISAZK3fpqjWwm2eiojoC9ylHpMCSJ6r9f6RAka30oVOkQLKKgvK3SIVjE/Y5DlA5Bsj3X9ygdgkWU1zT/HngnhzZKhyCZEFpU1164Tcfira91OGaDiIiIZGUzlQ0iIiJrIiCtK6T518LqMdkgIiKSgYDEbhSwG4WIiIjIKKxsEBERyUArbixStrcVTDaIiIhkwBlE67EbhYiIiGTFygYREZEMOINoPSYbREREMuAMovXYjUJERESyYmWDiIhIBpyuvB6TDSIiIhmwG6Uekw0iIiIZCHFjkbK9reCYDSIiIpIVKxtEREQy0EIFrYTnm0jZ1tow2SAiIpIBpyuvx24UIiIikhUrG0RERHKQOEDUlh6OwmSDiIhIBhyzUY/dKERERCQrVjaIiIhkwHk26jHZICIikgFnEK3HbhQiIiKSFSsbREREMuA8G/VsJtk4Vv05VKrmPXLXI7h5xw8AO7ZHKR2CRbxzYYfSIVjE2uAIpUOQbNePlUqHYBGfXf1M6RAkU7EYbhIBaXev2lCuYTvJBhERkTW5UdmQcOurDWUbTFOJiIhIVqxsEBERyYC3vtZjskFERCQD3vpaj90oREREJCtWNoiIiGTAbpR6TDaIiIhkwG6UeuxGISIiIlmxskFERCQDIXEGUXajEBERUZM4g2g9dqMQERHZkPT0dAQFBcHZ2RkRERHYv39/k+0//PBDBAcHw9nZGaGhodi+fbvutWvXrmHmzJkIDQ1FixYt4O/vjwkTJuDs2bMmxcRkg4iISAZ1D2KTsphq48aNSElJwdy5c1FQUICwsDBER0ejrKzMYPvc3FyMHTsWkyZNQmFhIWJiYhATE4MjR44AAKqrq1FQUIAXX3wRBQUF2Lx5M44dO4ZRo0aZFBeTDSIiIhnU3foqZTHV0qVLMXnyZCQmJqJ79+7IyMiAq6sr3n33XYPtly1bhmHDhmHGjBno1q0b5s+fj969e2P58uUAAA8PD2RnZyMuLg5du3bFXXfdheXLlyM/Px9FRUVGx8Vkg4iISAZaCywAUFlZqbfU1NQYPF5tbS3y8/MRFVX/9G07OztERUUhLy/P4DZ5eXl67QEgOjq60fYAUFFRAZVKhVatWjV5/n/GZIOIiMiKBQQEwMPDQ7ekpqYabFdeXg6NRgMfHx+99T4+PigpKTG4TUlJiUntr169ipkzZ2Ls2LFwd3c3+hx4NwoREZEMzB138eftAaC4uFjvD7uTk5PEyMxz7do1xMXFQQiBFStWmLQtkw0iIiIZWOrWV3d3d6OqCF5eXrC3t0dpaane+tLSUvj6+hrcxtfX16j2dYnGr7/+ii+//NKkqgbAbhQiIiKboFar0adPH+Tk5OjWabVa5OTkIDIy0uA2kZGReu0BIDs7W699XaJx4sQJfPHFF/D09DQ5NlY2iIiIZGCpbhRTpKSkICEhAX379kX//v2RlpaGqqoqJCYmAgAmTJiAdu3a6cZ9TJs2DYMHD8aSJUswYsQIZGVl4cCBA8jMzARwI9EYPXo0CgoK8Omnn0Kj0ejGc7Rp0wZqtdqouJhsEBERyUCJp77Gx8fj/PnzmDNnDkpKShAeHo4dO3boBoEWFRXBzq6+U2PAgAFYv349Zs+ejVmzZqFLly7YunUrQkJCAABnzpzBtm3bAADh4eF6x/rqq69w7733GhUXkw0iIiIbkpycjOTkZIOv7d69u8G62NhYxMbGGmwfFBQEYYGHtDDZICIikgEfMV+PyQYREZEMtJA4ZsNikShP8WQjNTUVmzdvxtGjR+Hi4oIBAwZg8eLF6Nq1q0n7cVH7QaVq3jfXLHxzmNIhSPbLZdv49dCK60qHYBGjFubcupGVu/KA4Wc6NDd2KuMG0lkznxZ9lQ5BMq24juraX5QO4y9H8b/Oe/bsQVJSEr799ltkZ2fj2rVrGDp0KKqqqpQOjYiIyGzCAoutULyysWPHDr1/r1mzBt7e3sjPz8egQYMUioqIiEgaIaR1hVhgXKbVUDzZuFlFRQWAG/fvGlJTU6P3EJrKysrbEhcREZEphJA4g6gNJRuKd6P8mVarxTPPPIOBAwfq7vG9WWpqqt4DaQICAm5zlERERGQKq0o2kpKScOTIEWRlZTXa5oUXXkBFRYVuKS4uvo0REhERGcdSj5i3BVbTjZKcnIxPP/0UX3/9Ndq3b99oOycnJ8WeeEdERGQsrQC0EjpSpNw2a20UTzaEEJgyZQq2bNmC3bt3o2PHjkqHRERERBakeLKRlJSE9evX4+OPP4abm5vuAS8eHh5wcXFRODoiIiLzWOoR87ZA8TEbK1asQEVFBe699174+fnplo0bNyodGhERkdnqnvoqZbEVilc2LPGAFyIiIrJeiicbREREtkj87z8p29sKJhtEREQy0EqcQdSWulEUH7NBREREto2VDSIiIhlInZiLk3oRERFRk4SQOGbDhm6gYLJBREQkA1Y26nHMBhEREcmKlQ0iIiIZsBulHpMNIiIiGQhI6wqxnVSD3ShEREQkM1Y2iIiIZKAVQuIj5m2ntsFkg4iISAacrrweu1GIiIhIVqxsEBERyYDzbNSzmWSjq30E7FVqpcOQ5MvzNUqHIFmlqlrpECyireOdSodgETXdWigdgmR9XVsqHYJFtIKL0iFIdrdn8/+TcVVbi0VVubflWFpIHLPBbhQiIiIi4zT/NJWIiMgK8W6Uekw2iIiIZMC7Ueox2SAiIpIBx2zU45gNIiIikhUrG0RERDJgZaMekw0iIiIZcMxGPXajEBERkaxY2SAiIpKBkNiNYkuVDSYbREREMtCqtFCpzJ90XGtDE5azG4WIiIhkxcoGERGRDLQQUPFuFABMNoiIiGQh/nfzq5TtbQW7UYiIiEhWrGwQERHJQAtI7EaxHaxsEBERyUCr0kpezJGeno6goCA4OzsjIiIC+/fvb7L9hx9+iODgYDg7OyM0NBTbt2/Xe33z5s0YOnQoPD09oVKpcPDgQZNjYrJBREQkA60F/jPVxo0bkZKSgrlz56KgoABhYWGIjo5GWVmZwfa5ubkYO3YsJk2ahMLCQsTExCAmJgZHjhzRtamqqsLdd9+NxYsXm30tmGwQERHZiKVLl2Ly5MlITExE9+7dkZGRAVdXV7z77rsG2y9btgzDhg3DjBkz0K1bN8yfPx+9e/fG8uXLdW0ee+wxzJkzB1FRUWbHxWSDiIhIBpaqbFRWVuotNTU1Bo9XW1uL/Px8vaTAzs4OUVFRyMvLM7hNXl5egyQiOjq60fbmYrJBREQkg7pbX6UsABAQEAAPDw/dkpqaavB45eXl0Gg08PHx0Vvv4+ODkpISg9uUlJSY1N5cvBuFiIjIihUXF8Pd3V33bycnJwWjMQ+TDSIiIhlY6tko7u7ueslGY7y8vGBvb4/S0lK99aWlpfD19TW4ja+vr0ntzcVuFCIiIhkIieM1TJ1BVK1Wo0+fPsjJydGt02q1yMnJQWRkpMFtIiMj9doDQHZ2dqPtzWUzlY0RbVvC2U6tdBiSZJQfVToEyezhqHQIFnHw8WNKh2ARW0beq3QIkp1WHVY6BIv4d2CY0iFI9sbpaqVDkEwjapUOQVYpKSlISEhA37590b9/f6SlpaGqqgqJiYkAgAkTJqBdu3a6cR/Tpk3D4MGDsWTJEowYMQJZWVk4cOAAMjMzdfv8448/UFRUhLNnzwIAjh278fno6+trdAXEZpINIiIiayKggZDQgSCgMXmb+Ph4nD9/HnPmzEFJSQnCw8OxY8cO3SDQoqIi2NnVxzRgwACsX78es2fPxqxZs9ClSxds3boVISEhujbbtm3TJSsAMGbMGADA3LlzMW/ePKPiYrJBREQkgxtjLqSP2TBVcnIykpOTDb62e/fuButiY2MRGxvb6P4mTpyIiRMnmhVLHY7ZICIiIlmxskFERCQDLQSkVTbMf4ibtWGyQUREJIMbYzZUkra3FUw2iIiIZKDUmA1rxDEbREREJCtWNoiIiGQgzJiY6+btbQWTDSIiIhlooQEkjNnQ2tCYDXajEBERkaxY2SAiIpIBu1HqMdkgIiKSgVZI7EYR7EYhIiIiMgorG0RERDJgN0o9JhtEREQyuJFsmN8VYkvJBrtRiIiISFasbBAREclACC20Up6NImynssFkg4iISAY3ukGkPIiNyQYRERE1QUi8dVXq9taEYzaIiIhIVqxsEBERyeDGiA12owBMNoiIiGRxY4AnB4gC7EYhIiIimbGyQUREJAMpE3pZYntrYjPJxitF70JKucoauDl3UToEya5e/0PpECzi8TWxSodgEZ9UZSkdgmRuTu2VDsEinji+XekQJEtq2/x/L2q0Nfj+yu05lhACkDJduRCWC0Zh7EYhIiIiWdlMZYOIiMiaSL2bhHejEBERUZNuTMplflcI70YhIiIiMhIrG0RERDKQWpmwpcoGkw0iIiIZcMxGPSYbREREMmBlo57VjNlIT09HUFAQnJ2dERERgf379ysdEhEREVmAVSQbGzduREpKCubOnYuCggKEhYUhOjoaZWVlSodGRERkFgGt5MVWWEWysXTpUkyePBmJiYno3r07MjIy4OrqinfffVfp0IiIiMwihEbyYisUTzZqa2uRn5+PqKgo3To7OztERUUhLy+vQfuamhpUVlbqLURERGS9FE82ysvLodFo4OPjo7fex8cHJSUlDdqnpqbCw8NDtwQEBNyuUImIiExQ92wUcxc+G0UxL7zwAioqKnRLcXGx0iERERE1IIRW8mIrFL/11cvLC/b29igtLdVbX1paCl9f3wbtnZyc4OTkdLvCIyIiIokUr2yo1Wr06dMHOTk5unVarRY5OTmIjIxUMDIiIiLz8W6UeopXNgAgJSUFCQkJ6Nu3L/r374+0tDRUVVUhMTFR6dCIiIjMpAWgkrC97YzZsIpkIz4+HufPn8ecOXNQUlKC8PBw7Nixo8GgUSIiImp+rCLZAIDk5GQkJycrHQYREZFlCImVDcHKBhERETVBSOxGEexGISIioqZxzEYdxe9GISIiItvGygYREZEshMTihO1UNphsEBERyULqqAsmG1ZD6EbrNv8fii084c9Wpte9JmqVDsEibOHnYQu/F4Bt/CxqtDVKhyBZrfbG77a4bXd6NP+/TZagErfvisvit99+48PYiIjIJMXFxWjfvr0s+7569So6duxo8GGipvL19cWpU6fg7OxsgciU0+yTDa1Wi7Nnz8LNzQ0qlZRRv42rrKxEQEAAiouL4e7uLssxbgdbOA9bOAfANs7DFs4B4HlYk9txDkIIXLp0Cf7+/rCzk+8eiatXr6K2VnqFVK1WN/tEA7CBbhQ7OzvZstObubu7N9tf4j+zhfOwhXMAbOM8bOEcAJ6HNZH7HDw8PGTbdx1nZ2ebSBIshbe+EhERkayYbBAREZGsmGwYwcnJCXPnzoWTk5PSoUhiC+dhC+cA2MZ52MI5ADwPa2IL50CGNfsBokRERGTdWNkgIiIiWTHZICIiIlkx2SAiIiJZMdkgIiIiWTHZuIX09HQEBQXB2dkZERER2L9/v9Ihmezrr7/GyJEj4e/vD5VKha1btyodkslSU1PRr18/uLm5wdvbGzExMTh27JjSYZlkxYoV6Nmzp27CosjISHz++edKhyXZokWLoFKp8MwzzygdiknmzZsHlUqltwQHBysdlsnOnDmD8ePHw9PTEy4uLggNDcWBAweUDsskQUFBDX4WKpUKSUlJSodGFsJkowkbN25ESkoK5s6di4KCAoSFhSE6OhplZWVKh2aSqqoqhIWFIT09XelQzLZnzx4kJSXh22+/RXZ2Nq5du4ahQ4eiqqpK6dCM1r59eyxatAj5+fk4cOAA7r//fjz00EP44YcflA7NbN999x1WrlyJnj17Kh2KWXr06IFz587plr179yodkkkuXLiAgQMHwtHREZ9//jl+/PFHLFmyBK1bt1Y6NJN89913ej+H7OxsAEBsbKzCkZHFCGpU//79RVJSku7fGo1G+Pv7i9TUVAWjkgaA2LJli9JhSFZWViYAiD179igdiiStW7cWb7/9ttJhmOXSpUuiS5cuIjs7WwwePFhMmzZN6ZBMMnfuXBEWFqZ0GJLMnDlT3H333UqHYXHTpk0TnTt3FlqtVulQyEJY2WhEbW0t8vPzERUVpVtnZ2eHqKgo5OXlKRgZAUBFRQUAoE2bNgpHYh6NRoOsrCxUVVUhMjJS6XDMkpSUhBEjRuj9jjQ3J06cgL+/Pzp16oRx48ahqKhI6ZBMsm3bNvTt2xexsbHw9vZGr169sGrVKqXDkqS2thbvv/8+Hn/8cdkerkm3H5ONRpSXl0Oj0cDHx0dvvY+Pj0UeG0zm02q1eOaZZzBw4ECEhIQoHY5JDh8+jJYtW8LJyQlPPfUUtmzZgu7duysdlsmysrJQUFCA1NRUpUMxW0REBNasWYMdO3ZgxYoVOHXqFO655x5cunRJ6dCM9ssvv2DFihXo0qULdu7ciaeffhpTp07F2rVrlQ7NbFu3bsXFixcxceJEpUMhC2r2T32lv56kpCQcOXKk2fWvA0DXrl1x8OBBVFRU4KOPPkJCQgL27NnTrBKO4uJiTJs2DdnZ2c36qZbDhw/X/X/Pnj0RERGBwMBAfPDBB5g0aZKCkRlPq9Wib9++WLhwIQCgV69eOHLkCDIyMpCQkKBwdOZ55513MHz4cPj7+ysdClkQKxuN8PLygr29PUpLS/XWl5aWwtfXV6GoKDk5GZ9++im++uortG/fXulwTKZWq3HHHXegT58+SE1NRVhYGJYtW6Z0WCbJz89HWVkZevfuDQcHBzg4OGDPnj1444034ODgAI1Go3SIZmnVqhXuvPNOnDx5UulQjObn59cgUe3WrVuz6w6q8+uvv+KLL77AP/7xD6VDIQtjstEItVqNPn36ICcnR7dOq9UiJyen2faxN2dCCCQnJ2PLli348ssv0bFjR6VDsgitVouamhqlwzDJkCFDcPjwYRw8eFC39O3bF+PGjcPBgwdhb2+vdIhmuXz5Mn7++Wf4+fkpHYrRBg4c2OAW8OPHjyMwMFChiKRZvXo1vL29MWLECKVDIQtjN0oTUlJSkJCQgL59+6J///5IS0tDVVUVEhMTlQ7NJJcvX9b7tnbq1CkcPHgQbdq0QYcOHRSMzHhJSUlYv349Pv74Y7i5uenGzXh4eMDFxUXh6IzzwgsvYPjw4ejQoQMuXbqE9evXY/fu3di5c6fSoZnEzc2twViZFi1awNPTs1mNoZk+fTpGjhyJwMBAnD17FnPnzoW9vT3Gjh2rdGhGe/bZZzFgwAAsXLgQcXFx2L9/PzIzM5GZmal0aCbTarVYvXo1EhIS4ODAP002R+nbYazdm2++KTp06CDUarXo37+/+Pbbb5UOyWRfffWVANBgSUhIUDo0oxmKH4BYvXq10qEZ7fHHHxeBgYFCrVaLtm3biiFDhohdu3YpHZZFNMdbX+Pj44Wfn59Qq9WiXbt2Ij4+Xpw8eVLpsEz2ySefiJCQEOHk5CSCg4NFZmam0iGZZefOnQKAOHbsmNKhkAz4iHkiIiKSFcdsEBERkayYbBAREZGsmGwQERGRrJhsEBERkayYbBAREZGsmGwQERGRrJhsEBERkayYbBAREZGsmGwQERGRrJhsEDVT9957L5555hlZ9j1jxgyMHDlSln0T0V8Ppysnaqb++OMPODo6ws3NDcCN5CM8PBxpaWmS933x4kXY29vr9k1EJAUfrUfUTLVp00a2fbdq1Uq2fRPRXw+7UYis2EcffYTQ0FC4uLjA09MTUVFRqKqqAqDfjTJx4kTs2bMHy5Ytg0qlgkqlwunTp6HVapGamoqOHTvCxcUFYWFh+Oijj5o8Znl5OVQqFY4cOSL36RHRXwQrG0RW6ty5cxg7dixeffVVPPzww7h06RK++eYbGOr5XLZsGY4fP46QkBC8/PLLAIC2bdsiNTUV77//PjIyMtClSxd8/fXXGD9+PNq2bYvBgwcbPO73338PJycnBAcHy3p+RPTXwWSDyEqdO3cO169fxyOPPILAwEAAQGhoqMG2Hh4eUKvVcHV1ha+vLwCgpqYGCxcuxBdffIHIyEgAQKdOnbB3716sXLmyyWSjR48ecHDgxwMRWQY/TYisVFhYGIYMGYLQ0FBER0dj6NChGD16NFq3bm3U9idPnkR1dTX+9re/6a2vra1Fr169Gt3u4MGDCA8PlxI6EZEeJhtEVsre3h7Z2dnIzc3Frl278Oabb+Jf//oX9u3bh44dO95y+8uXLwMAPvvsM7Rr107vNScnp0a3+/777zFp0iRpwRMR/QkHiBJZMZVKhYEDB+Kll15CYWEh1Go1tmzZYrCtWq2GRqPR/bt79+5wcnJCUVER7rjjDr0lICDA4D5qa2vx008/ISwsTJbzIaK/JlY2iKzUvn37kJOTg6FDh8Lb2xv79u3D+fPn0a1bN4Ptg4KCsG/fPpw+fRotW7ZEmzZtMH36dDz77LPQarW4++67UVFRgf/+979wd3dHQkJCg3389NNPuHbtGpMNIrIoJhtEVsrd3R1ff/010tLSUFlZicDAQCxZsgTDhw832H769OlISEhA9+7dceXKFZw6dQrz58/X3ZXyyy+/oFWrVujduzdmzZplcB8HDx5EYGAg59kgIoviDKJEpJOcnIyysjJ88MEHSodCRDaEYzaICFevXkV+fj42bdqE6OhopcMhIhvDZIOIkJaWhqioKDz00EOYMGGC0uEQkY1hNwoRERHJipUNIiIikhWTDSIiIpIVkw0iIiKSFZMNIiIikhWTDSIiIpIVkw0iIiKSFZMNIiIikhWTDSIiIpIVkw0iIiKSFZMNIiIikhWTDSIiIpLV/wMhfHsEXlmRxQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0, N_s-1, N_s)\n",
    "y = np.linspace(0, (N_t - 1) * dt, N_t)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "perc_errors = []\n",
    "for t in range(N_t):\n",
    "    perc_error_t = []\n",
    "    for j in range(N_s):\n",
    "        perc_error = np.abs((qc_occ_nums[t][j] - occs[t][j]))\n",
    "        perc_error_t.append(perc_error)\n",
    "    perc_errors.append(perc_error_t)\n",
    "\n",
    "plt.pcolormesh(X, Y, perc_errors, cmap='inferno')\n",
    "plt.xlabel(r'site $j$')\n",
    "plt.ylabel(r'time $t$')\n",
    "title = r'(c) Absolute Difference: $N_s = $' + str(N_s)\n",
    "title += r', $J = $' + str(J)\n",
    "# title += r', $h = $' + str(h)\n",
    "title += r', $g = $' + str(g)\n",
    "title += r', $\\Delta t=$' + str(dt)\n",
    "plt.title(title)\n",
    "plt.colorbar()\n",
    "# plt.savefig('qc-occ-nums.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.06500510701066574"
      ]
     },
     "execution_count": 128,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(perc_errors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "nks = np.arange(-N_s + 1, N_s + 1, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "ETADAG8 = []\n",
    "for i in range(N_s):\n",
    "    eta = np.zeros((2**N_s, 2**N_s), dtype=complex)\n",
    "    for j in range(N_s):\n",
    "        arg = 1 / np.sqrt(N_s) * np.exp(-1j * nks[i] * np.pi / N_s * j)\n",
    "        eta += arg * CDAG[j]\n",
    "    eta /= np.linalg.norm(eta)\n",
    "    ETADAG8.append(eta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "ETADAG12 = []\n",
    "for i in range(N_s):\n",
    "    eta = ETADAG8[i]\n",
    "    eta = np.kron(np.identity(2**4, dtype=complex), eta)\n",
    "    eta /= np.linalg.norm(eta)\n",
    "    ETADAG12.append(eta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.quantum_info import Statevector\n",
    "qreg = QuantumRegister(N_s+4)\n",
    "circ = QuantumCircuit(qreg)\n",
    "\n",
    "# State-Prep\n",
    "## Ground state preparation\n",
    "circ.prepare_state(ground_state, range(N_s)) # Exact\n",
    "ground_state12 = Statevector.from_instruction(circ).reverse_qargs().data\n",
    "\n",
    "circ.x(N_s) # Prepare control\n",
    "circ.x(N_s+2) # Prepare control\n",
    "\n",
    "## Right moving wave packet\n",
    "### V^dag\n",
    "for i in range(N_s//2):\n",
    "    circ.rz(betas_r[i], i)\n",
    "for i in range(N_s//2-1):\n",
    "    j = N_s//2 - 1 - i\n",
    "    circ.append(givens_circ(angles_r[i]), [j-1, j])\n",
    "###\n",
    "circ.ccx(N_s, 0, N_s+1) # Remove |1> part of system qubit\n",
    "circ.x(0) # Excite system qubit\n",
    "### V\n",
    "for i in range(N_s//2-1):\n",
    "    j = N_s//2 - 2 - i\n",
    "    circ.append(givens_circ(-angles_r[j]), [i, i+1])\n",
    "for i in range(N_s//2):\n",
    "    circ.rz(-betas_r[i], i)\n",
    "\n",
    "# # Left Moving Wave Packet\n",
    "# # V^dag\n",
    "for i in range(N_s//2):\n",
    "    circ.rz(betas_l[i], i+N_s//2)\n",
    "for i in range(N_s//2-1):\n",
    "    j = N_s - 1 - i\n",
    "    circ.append(givens_circ(angles_l[i]), [j-1, j])\n",
    "###\n",
    "circ.ccx(N_s+2, N_s//2, N_s+3) # Remove |1> part of system qubit\n",
    "for i in range(N_s//2): # Obey anticommutation relations\n",
    "    circ.z(i)\n",
    "circ.x(N_s//2) # Excite system qubit\n",
    "### V\n",
    "for i in range(N_s//2-1):\n",
    "    j = N_s//2 - 2 - i\n",
    "    circ.append(givens_circ(-angles_l[j]), [i+N_s//2, i+N_s//2+1])\n",
    "for i in range(N_s//2):\n",
    "    circ.rz(-betas_l[i], i+N_s//2)\n",
    "\n",
    "state = Statevector.from_instruction(circ).reverse_qargs().data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "ground_state_12 = ground_state12.reshape((2**12))\n",
    "state = state.reshape((2**12))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "four_spec_8 = []\n",
    "four_spec_12 = []\n",
    "for i in range(N_s-1):\n",
    "    for j in range(i+1, N_s):\n",
    "        k_state8 = ETADAG8[i] @ ETADAG8[j] @ ground_state\n",
    "        k_state8 /= np.linalg.norm(k_state8)\n",
    "        # k_state12 = ETADAG12[i] @ ETADAG12[j] @ ground_state_12\n",
    "        # k_state12 /= np.linalg.norm(k_state12)\n",
    "        \n",
    "        amp8 = np.inner(k_state8.conj(), exact_init_state)\n",
    "        # amp12 = np.inner(k_state12.conj(), state)\n",
    "\n",
    "        four_spec_8.append(amp8 * amp8.conj())\n",
    "        # four_spec_12.append(amp12 * amp12.conj())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "four_spec_12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/site-packages/matplotlib/cbook.py:1762: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return math.isfinite(val)\n",
      "/Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/site-packages/matplotlib/transforms.py:767: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  points = np.asarray(points, float)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<BarContainer object of 0 artists>"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAo5klEQVR4nO3df1DU953H8RegLBqFRDlZQeJqtCpVQUEoJhdyk51Aa5uQeIR4XqXUIZOeJJq9oxFP4Xr2ujZRBxuZcHbGtJ0LlXOuWpt45LyNmOuJ4YQwOWPOpF5SOMku2F6gYgMO+70/cllnK6CLP/YDPh8z3wl89v357vv7zTf6yofvd4mwLMsSAACAwSLD3QAAAMDVEFgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYbF+4GbhS/36+Ojg5NnjxZERER4W4HAABcA8uy9Lvf/U6JiYmKjBx6HWXMBJaOjg4lJyeHuw0AADAC7e3tmjFjxpCvj5nAMnnyZEmfHXBsbGyYuwEAANeip6dHycnJgb/HhzJmAsvnPwaKjY0lsAAAMMpc7XYObroFAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMN64cDcAALeKY+NrIc/5aNuKm9AJgFCxwgIAAIw3osBSXV0th8OhmJgYZWVlqampacjad999VytXrpTD4VBERISqqqoGrTt37pz+/M//XFOnTtWECRO0aNEinTx5ciTtAQCAMSbkwFJXVyeXy6XKykq1tLQoNTVVubm56uzsHLT+4sWLmj17trZt2ya73T5ozf/+7//q3nvv1fjx4/XP//zPOn36tHbs2KG77ror1PYAAMAYFPI9LDt37lRJSYmKi4slSTU1NXrttde0d+9ebdy48Yr6ZcuWadmyZZI06OuS9P3vf1/Jycl6+eWXA2OzZs0KtTUAADBGhbTC0t/fr+bmZjmdzss7iIyU0+lUY2PjiJs4dOiQMjIyVFBQoGnTpmnJkiX64Q9/OOycvr4+9fT0BG0AAGBsCimwnD9/XgMDA0pISAgaT0hIkNfrHXET//3f/62XXnpJc+fO1euvv65vfetbeuaZZ/TjH/94yDlut1txcXGBLTk5ecTvDwAAzGbEU0J+v19Lly7V9773PS1ZskRPPvmkSkpKVFNTM+Sc8vJydXd3B7b29vZb2DEAALiVQgos8fHxioqKks/nCxr3+XxD3lB7LaZPn66UlJSgsQULFqitrW3IOTabTbGxsUEbAAAYm0IKLNHR0UpPT5fH4wmM+f1+eTweZWdnj7iJe++9V2fOnAkae//99zVz5swR7xMAAIwdIT8l5HK5VFRUpIyMDGVmZqqqqkq9vb2Bp4bWrFmjpKQkud1uSZ/dqHv69OnA1+fOnVNra6smTZqkOXPmSJKeffZZLV++XN/73vf0+OOPq6mpSXv27NGePXtu1HECAIBRLOTAUlhYqK6uLlVUVMjr9SotLU319fWBG3Hb2toUGXl54aajo0NLliwJfL99+3Zt375dOTk5amhokPTZo88HDhxQeXm5/vZv/1azZs1SVVWVVq9efZ2HBwAAxoIIy7KscDdxI/T09CguLk7d3d3czwJgUPwuIcA81/r3txFPCQEAAAyHwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMN6IAkt1dbUcDodiYmKUlZWlpqamIWvfffddrVy5Ug6HQxEREaqqqhp239u2bVNERIQ2bNgwktYAAMAYFHJgqaurk8vlUmVlpVpaWpSamqrc3Fx1dnYOWn/x4kXNnj1b27Ztk91uH3bf//Ef/6G///u/1+LFi0NtCwAAjGEhB5adO3eqpKRExcXFSklJUU1NjSZOnKi9e/cOWr9s2TK98MILeuKJJ2Sz2Ybc74ULF7R69Wr98Ic/1F133RVqWwAAYAwLKbD09/erublZTqfz8g4iI+V0OtXY2Hhdjaxbt04rVqwI2vdw+vr61NPTE7QBAICxKaTAcv78eQ0MDCghISFoPCEhQV6vd8RN7Nu3Ty0tLXK73dc8x+12Ky4uLrAlJyeP+P0BAIDZwv6UUHt7u9avX69XXnlFMTEx1zyvvLxc3d3dga29vf0mdgkAAMJpXCjF8fHxioqKks/nCxr3+XxXvaF2KM3Nzers7NTSpUsDYwMDA3rzzTe1e/du9fX1KSoq6op5Nptt2HtiAADA2BHSCkt0dLTS09Pl8XgCY36/Xx6PR9nZ2SNq4MEHH9R//ud/qrW1NbBlZGRo9erVam1tHTSsAACA20tIKyyS5HK5VFRUpIyMDGVmZqqqqkq9vb0qLi6WJK1Zs0ZJSUmB+1H6+/t1+vTpwNfnzp1Ta2urJk2apDlz5mjy5MlauHBh0Hvccccdmjp16hXjAADg9hRyYCksLFRXV5cqKirk9XqVlpam+vr6wI24bW1tioy8vHDT0dGhJUuWBL7fvn27tm/frpycHDU0NFz/EQAAgDEvwrIsK9xN3Ag9PT2Ki4tTd3e3YmNjw90OAAM5Nr4W8pyPtq24CZ0A+Ny1/v0d9qeEAAAArobAAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYb0SBpbq6Wg6HQzExMcrKylJTU9OQte+++65Wrlwph8OhiIgIVVVVXVHjdru1bNkyTZ48WdOmTVN+fr7OnDkzktYAAMAYFHJgqaurk8vlUmVlpVpaWpSamqrc3Fx1dnYOWn/x4kXNnj1b27Ztk91uH7Tm2LFjWrdunU6cOKEjR47o0qVLeuihh9Tb2xtqewAAYAyKsCzLCmVCVlaWli1bpt27d0uS/H6/kpOT9fTTT2vjxo3DznU4HNqwYYM2bNgwbF1XV5emTZumY8eO6f7777+mvnp6ehQXF6fu7m7FxsZe0xwAtxfHxtdCnvPRthU3oRMAn7vWv79DWmHp7+9Xc3OznE7n5R1ERsrpdKqxsXHk3f6B7u5uSdKUKVOGrOnr61NPT0/QBgAAxqaQAsv58+c1MDCghISEoPGEhAR5vd4b0pDf79eGDRt07733auHChUPWud1uxcXFBbbk5OQb8v4AAMA8xj0ltG7dOp06dUr79u0btq68vFzd3d2Brb29/RZ1CAAAbrVxoRTHx8crKipKPp8vaNzn8w15Q20oSktL9eqrr+rNN9/UjBkzhq212Wyy2WzX/Z4AAMB8Ia2wREdHKz09XR6PJzDm9/vl8XiUnZ094iYsy1JpaakOHDigN954Q7NmzRrxvgAAwNgT0gqLJLlcLhUVFSkjI0OZmZmqqqpSb2+viouLJUlr1qxRUlKS3G63pM9u1D19+nTg63Pnzqm1tVWTJk3SnDlzJH32Y6Da2lr9/Oc/1+TJkwP3w8TFxWnChAk35EABAMDoFXJgKSwsVFdXlyoqKuT1epWWlqb6+vrAjbhtbW2KjLy8cNPR0aElS5YEvt++fbu2b9+unJwcNTQ0SJJeeuklSdIDDzwQ9F4vv/yyvvGNb4TaIgAAGGNC/hwWU/E5LACuhs9hAcxzUz6HBQAAIBwILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA440osFRXV8vhcCgmJkZZWVlqamoasvbdd9/VypUr5XA4FBERoaqqquveJwAAuL2EHFjq6urkcrlUWVmplpYWpaamKjc3V52dnYPWX7x4UbNnz9a2bdtkt9tvyD4BAMDtJeTAsnPnTpWUlKi4uFgpKSmqqanRxIkTtXfv3kHrly1bphdeeEFPPPGEbDbbDdknAAC4vYQUWPr7+9Xc3Cyn03l5B5GRcjqdamxsHFEDI91nX1+fenp6gjYAADA2hRRYzp8/r4GBASUkJASNJyQkyOv1jqiBke7T7XYrLi4usCUnJ4/o/QEAgPlG7VNC5eXl6u7uDmzt7e3hbgkAANwk40Ipjo+PV1RUlHw+X9C4z+cb8obam7VPm8025D0xwFjk2PhayHM+2rbiJnQSPiacAxN6uB6jvX/cvkJaYYmOjlZ6ero8Hk9gzO/3y+PxKDs7e0QN3Ix9AgCAsSWkFRZJcrlcKioqUkZGhjIzM1VVVaXe3l4VFxdLktasWaOkpCS53W5Jn91Ue/r06cDX586dU2trqyZNmqQ5c+Zc0z4BAMDtLeTAUlhYqK6uLlVUVMjr9SotLU319fWBm2bb2toUGXl54aajo0NLliwJfL99+3Zt375dOTk5amhouKZ9AgCA21vIgUWSSktLVVpaOuhrn4eQzzkcDlmWdV37BAAAt7dR+5QQAAC4fRBYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOONC3cDAIBr59j4WshzPtq24iZ0AtxarLAAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGG9EgaW6uloOh0MxMTHKyspSU1PTsPX79+/X/PnzFRMTo0WLFunw4cNBr1+4cEGlpaWaMWOGJkyYoJSUFNXU1IykNQAAMAaFHFjq6urkcrlUWVmplpYWpaamKjc3V52dnYPWHz9+XKtWrdLatWv19ttvKz8/X/n5+Tp16lSgxuVyqb6+Xv/wD/+g9957Txs2bFBpaakOHTo08iMDAABjRsiBZefOnSopKVFxcXFgJWTixInau3fvoPW7du1SXl6eysrKtGDBAm3dulVLly7V7t27AzXHjx9XUVGRHnjgATkcDj355JNKTU296soNAAC4PYQUWPr7+9Xc3Cyn03l5B5GRcjqdamxsHHROY2NjUL0k5ebmBtUvX75chw4d0rlz52RZlo4ePar3339fDz300JC99PX1qaenJ2gDAABjU0iB5fz58xoYGFBCQkLQeEJCgrxe76BzvF7vVetffPFFpaSkaMaMGYqOjlZeXp6qq6t1//33D9mL2+1WXFxcYEtOTg7lUAAAwChixFNCL774ok6cOKFDhw6publZO3bs0Lp16/Sv//qvQ84pLy9Xd3d3YGtvb7+FHQMAgFtpXCjF8fHxioqKks/nCxr3+Xyy2+2DzrHb7cPW//73v9emTZt04MABrVixQpK0ePFitba2avv27Vf8OOlzNptNNpstlPYBAMAoFdIKS3R0tNLT0+XxeAJjfr9fHo9H2dnZg87Jzs4OqpekI0eOBOovXbqkS5cuKTIyuJWoqCj5/f5Q2gMAAGNUSCss0mePIBcVFSkjI0OZmZmqqqpSb2+viouLJUlr1qxRUlKS3G63JGn9+vXKycnRjh07tGLFCu3bt08nT57Unj17JEmxsbHKyclRWVmZJkyYoJkzZ+rYsWP6yU9+op07d97AQwUAAKNVyIGlsLBQXV1dqqiokNfrVVpamurr6wM31ra1tQWtlixfvly1tbXavHmzNm3apLlz5+rgwYNauHBhoGbfvn0qLy/X6tWr9dvf/lYzZ87U3/3d3+mpp566AYcIAABGu5ADiySVlpaqtLR00NcaGhquGCsoKFBBQcGQ+7Pb7Xr55ZdH0goAALgNGPGUEAAAwHAILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA440osFRXV8vhcCgmJkZZWVlqamoatn7//v2aP3++YmJitGjRIh0+fPiKmvfee08PP/yw4uLidMcdd2jZsmVqa2sbSXsAAGCMCTmw1NXVyeVyqbKyUi0tLUpNTVVubq46OzsHrT9+/LhWrVqltWvX6u2331Z+fr7y8/N16tSpQM3Zs2d13333af78+WpoaNA777yjLVu2KCYmZuRHBgAAxoyQA8vOnTtVUlKi4uJipaSkqKamRhMnTtTevXsHrd+1a5fy8vJUVlamBQsWaOvWrVq6dKl2794dqPnrv/5rfeUrX9Hzzz+vJUuW6J577tHDDz+sadOmjfzIAADAmBFSYOnv71dzc7OcTuflHURGyul0qrGxcdA5jY2NQfWSlJubG6j3+/167bXX9IUvfEG5ubmaNm2asrKydPDgwWF76evrU09PT9AGAADGppACy/nz5zUwMKCEhISg8YSEBHm93kHneL3eYes7Ozt14cIFbdu2TXl5efqXf/kXPfroo3rsscd07NixIXtxu92Ki4sLbMnJyaEcCgAAGEXC/pSQ3++XJD3yyCN69tlnlZaWpo0bN+qrX/2qampqhpxXXl6u7u7uwNbe3n6rWgYAALfYuFCK4+PjFRUVJZ/PFzTu8/lkt9sHnWO324etj4+P17hx45SSkhJUs2DBAv3yl78cshebzSabzRZK+wAAYJQKaYUlOjpa6enp8ng8gTG/3y+Px6Ps7OxB52RnZwfVS9KRI0cC9dHR0Vq2bJnOnDkTVPP+++9r5syZobQHAADGqJBWWCTJ5XKpqKhIGRkZyszMVFVVlXp7e1VcXCxJWrNmjZKSkuR2uyVJ69evV05Ojnbs2KEVK1Zo3759OnnypPbs2RPYZ1lZmQoLC3X//ffrT/7kT1RfX69f/OIXamhouDFHCRjAsfG1kOd8tG3FTegEuD5cywiHkANLYWGhurq6VFFRIa/Xq7S0NNXX1wdurG1ra1Nk5OWFm+XLl6u2tlabN2/Wpk2bNHfuXB08eFALFy4M1Dz66KOqqamR2+3WM888o3nz5umf/umfdN99992AQwQAAKNdyIFFkkpLS1VaWjroa4OtihQUFKigoGDYfX7zm9/UN7/5zZG0A+Aa8X/GAEarsD8lBAAAcDUEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeCMKLNXV1XI4HIqJiVFWVpaampqGrd+/f7/mz5+vmJgYLVq0SIcPHx6y9qmnnlJERISqqqpG0hoAABiDQg4sdXV1crlcqqysVEtLi1JTU5Wbm6vOzs5B648fP65Vq1Zp7dq1evvtt5Wfn6/8/HydOnXqitoDBw7oxIkTSkxMDP1IAADAmBVyYNm5c6dKSkpUXFyslJQU1dTUaOLEidq7d++g9bt27VJeXp7Kysq0YMECbd26VUuXLtXu3buD6s6dO6enn35ar7zyisaPHz+yowEAAGNSSIGlv79fzc3Ncjqdl3cQGSmn06nGxsZB5zQ2NgbVS1Jubm5Qvd/v19e//nWVlZXpi1/84jX10tfXp56enqANAACMTSEFlvPnz2tgYEAJCQlB4wkJCfJ6vYPO8Xq9V63//ve/r3HjxumZZ5655l7cbrfi4uICW3JycghHAgAARpOwPyXU3NysXbt26Uc/+pEiIiKueV55ebm6u7sDW3t7+03sEgAAhFNIgSU+Pl5RUVHy+XxB4z6fT3a7fdA5drt92Pp/+7d/U2dnp+6++26NGzdO48aN069//Wv95V/+pRwOx5C92Gw2xcbGBm0AAGBsCimwREdHKz09XR6PJzDm9/vl8XiUnZ096Jzs7Oygekk6cuRIoP7rX/+63nnnHbW2tga2xMRElZWV6fXXXw/1eAAAwBg0LtQJLpdLRUVFysjIUGZmpqqqqtTb26vi4mJJ0po1a5SUlCS32y1JWr9+vXJycrRjxw6tWLFC+/bt08mTJ7Vnzx5J0tSpUzV16tSg9xg/frzsdrvmzZt3vccHAADGgJADS2Fhobq6ulRRUSGv16u0tDTV19cHbqxta2tTZOTlhZvly5ertrZWmzdv1qZNmzR37lwdPHhQCxcuvHFHAQAAxrSQA4sklZaWqrS0dNDXGhoarhgrKChQQUHBNe//o48+GklbAABgjAr7U0IAAABXQ2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGC8ceFuAMDtw7HxtZDnfLRtxU3oBKMZ19HtiRUWAABgvBEFlurqajkcDsXExCgrK0tNTU3D1u/fv1/z589XTEyMFi1apMOHDwdeu3Tpkp577jktWrRId9xxhxITE7VmzRp1dHSMpDUAADAGhRxY6urq5HK5VFlZqZaWFqWmpio3N1ednZ2D1h8/flyrVq3S2rVr9fbbbys/P1/5+fk6deqUJOnixYtqaWnRli1b1NLSop/97Gc6c+aMHn744es7MgAAMGaEHFh27typkpISFRcXKyUlRTU1NZo4caL27t07aP2uXbuUl5ensrIyLViwQFu3btXSpUu1e/duSVJcXJyOHDmixx9/XPPmzdOXvvQl7d69W83NzWpra7u+owMAAGNCSIGlv79fzc3Ncjqdl3cQGSmn06nGxsZB5zQ2NgbVS1Jubu6Q9ZLU3d2tiIgI3XnnnUPW9PX1qaenJ2gDAABjU0iB5fz58xoYGFBCQkLQeEJCgrxe76BzvF5vSPWffvqpnnvuOa1atUqxsbFD9uJ2uxUXFxfYkpOTQzkUAAAwihj1lNClS5f0+OOPy7IsvfTSS8PWlpeXq7u7O7C1t7ffoi4BAMCtFtLnsMTHxysqKko+ny9o3OfzyW63DzrHbrdfU/3nYeXXv/613njjjWFXVyTJZrPJZrOF0j4AABilQgos0dHRSk9Pl8fjUX5+viTJ7/fL4/GotLR00DnZ2dnyeDzasGFDYOzIkSPKzs4OfP95WPnggw909OhRTZ06NfQjwU3FBzUBNwb/LQEjE/In3bpcLhUVFSkjI0OZmZmqqqpSb2+viouLJUlr1qxRUlKS3G63JGn9+vXKycnRjh07tGLFCu3bt08nT57Unj17JH0WVv70T/9ULS0tevXVVzUwMBC4v2XKlCmKjo6+UccKAABGqZADS2Fhobq6ulRRUSGv16u0tDTV19cHbqxta2tTZOTlW2OWL1+u2tpabd68WZs2bdLcuXN18OBBLVy4UJJ07tw5HTp0SJKUlpYW9F5Hjx7VAw88MMJDAwAAY8WIfpdQaWnpkD8CamhouGKsoKBABQUFg9Y7HA5ZljWSNgAAwG3CqKeEAAAABkNgAQAAxiOwAAAA4xFYAACA8UZ00y1wO+LzM4Cxgf+WRydWWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjMdjzbgleIxwbODfI4BwYYUFAAAYjxWW2wT/ZwzAFPx5hJFghQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8fhtzaMEv90UAMzAn8fhwQoLAAAwHoEFAAAYj8ACAACMxz0sGBX4mTEA3N5GtMJSXV0th8OhmJgYZWVlqampadj6/fv3a/78+YqJidGiRYt0+PDhoNcty1JFRYWmT5+uCRMmyOl06oMPPhhJawAAYAwKeYWlrq5OLpdLNTU1ysrKUlVVlXJzc3XmzBlNmzbtivrjx49r1apVcrvd+upXv6ra2lrl5+erpaVFCxculCQ9//zz+sEPfqAf//jHmjVrlrZs2aLc3FydPn1aMTEx13+UYcbqAADgRrod/14JeYVl586dKikpUXFxsVJSUlRTU6OJEydq7969g9bv2rVLeXl5Kisr04IFC7R161YtXbpUu3fvlvTZ6kpVVZU2b96sRx55RIsXL9ZPfvITdXR06ODBg9d1cAAAYGwIaYWlv79fzc3NKi8vD4xFRkbK6XSqsbFx0DmNjY1yuVxBY7m5uYEw8uGHH8rr9crpdAZej4uLU1ZWlhobG/XEE08Mut++vj719fUFvu/u7pYk9fT0hHJI12Rh5eshzzn1ndzA1/6+iyHP/8PjuN593O7zpfD/ewz3fBN6GO3zTegh3PNN6GG0z5fC/+fR9b7/jfR5X5ZlDV9oheDcuXOWJOv48eNB42VlZVZmZuagc8aPH2/V1tYGjVVXV1vTpk2zLMuy/v3f/92SZHV0dATVFBQUWI8//viQvVRWVlqS2NjY2NjY2MbA1t7ePmwGGbVPCZWXlwet3Pj9fv32t7/V1KlTFRERcUt66OnpUXJystrb2xUbG3tL3nOs4RxeP87h9eMcXj/O4fW7Xc+hZVn63e9+p8TExGHrQgos8fHxioqKks/nCxr3+Xyy2+2DzrHb7cPWf/5Pn8+n6dOnB9WkpaUN2YvNZpPNZgsau/POO6/1UG6o2NjY2+riuhk4h9ePc3j9OIfXj3N4/W7HcxgXF3fVmpBuuo2OjlZ6ero8Hk9gzO/3y+PxKDs7e9A52dnZQfWSdOTIkUD9rFmzZLfbg2p6enr01ltvDblPAABwewn5R0Iul0tFRUXKyMhQZmamqqqq1Nvbq+LiYknSmjVrlJSUJLfbLUlav369cnJytGPHDq1YsUL79u3TyZMntWfPHklSRESENmzYoO9+97uaO3du4LHmxMRE5efn37gjBQAAo1bIgaWwsFBdXV2qqKiQ1+tVWlqa6uvrlZCQIElqa2tTZOTlhZvly5ertrZWmzdv1qZNmzR37lwdPHgw8BkskvTtb39bvb29evLJJ/XJJ5/ovvvuU319vfGfwWKz2VRZWXnFj6Zw7TiH149zeP04h9ePc3j9OIfDi7Csqz1HBAAAEF788kMAAGA8AgsAADAegQUAABiPwAIAAIxHYBmh6upqORwOxcTEKCsrS01NTeFuadT4m7/5G0VERARt8+fPD3dbxnvzzTf1ta99TYmJiYqIiLjil4NalqWKigpNnz5dEyZMkNPp1AcffBCeZg11tXP4jW9844prMy8vLzzNGsjtdmvZsmWaPHmypk2bpvz8fJ05cyao5tNPP9W6des0depUTZo0SStXrrziw0NvZ9dyDh944IErrsOnnnoqTB2bg8AyAnV1dXK5XKqsrFRLS4tSU1OVm5urzs7OcLc2anzxi1/Uxx9/HNh++ctfhrsl4/X29io1NVXV1dWDvv7888/rBz/4gWpqavTWW2/pjjvuUG5urj799NNb3Km5rnYOJSkvLy/o2vzpT396Czs027Fjx7Ru3TqdOHFCR44c0aVLl/TQQw+pt7c3UPPss8/qF7/4hfbv369jx46po6NDjz32WBi7Nsu1nENJKikpCboOn3/++TB1bJBhf9MQBpWZmWmtW7cu8P3AwICVmJhoud3uMHY1elRWVlqpqanhbmNUk2QdOHAg8L3f77fsdrv1wgsvBMY++eQTy2azWT/96U/D0KH5/vAcWpZlFRUVWY888khY+hmNOjs7LUnWsWPHLMv67JobP368tX///kDNe++9Z0myGhsbw9Wm0f7wHFqWZeXk5Fjr168PX1OGYoUlRP39/WpubpbT6QyMRUZGyul0qrGxMYydjS4ffPCBEhMTNXv2bK1evVptbW3hbmlU+/DDD+X1eoOuy7i4OGVlZXFdhqihoUHTpk3TvHnz9K1vfUu/+c1vwt2Ssbq7uyVJU6ZMkSQ1Nzfr0qVLQdfh/Pnzdffdd3MdDuEPz+HnXnnlFcXHx2vhwoUqLy/XxYsXw9GeUUbtb2sOl/Pnz2tgYCDwyb6fS0hI0H/913+FqavRJSsrSz/60Y80b948ffzxx/rOd76jP/7jP9apU6c0efLkcLc3Knm9Xkka9Lr8/DVcXV5enh577DHNmjVLZ8+e1aZNm/TlL39ZjY2NioqKCnd7RvH7/dqwYYPuvffewCeXe71eRUdHX/GLaLkOBzfYOZSkP/uzP9PMmTOVmJiod955R88995zOnDmjn/3sZ2HsNvwILLjlvvzlLwe+Xrx4sbKysjRz5kz94z/+o9auXRvGznC7e+KJJwJfL1q0SIsXL9Y999yjhoYGPfjgg2HszDzr1q3TqVOnuP/sOgx1Dp988snA14sWLdL06dP14IMP6uzZs7rnnntudZvG4EdCIYqPj1dUVNQVd737fD7Z7fYwdTW63XnnnfrCF76gX/3qV+FuZdT6/NrjuryxZs+erfj4eK7NP1BaWqpXX31VR48e1YwZMwLjdrtd/f39+uSTT4LquQ6vNNQ5HExWVpYk3fbXIYElRNHR0UpPT5fH4wmM+f1+eTweZWdnh7Gz0evChQs6e/aspk+fHu5WRq1Zs2bJbrcHXZc9PT166623uC6vw//8z//oN7/5Ddfm/7MsS6WlpTpw4IDeeOMNzZo1K+j19PR0jR8/Pug6PHPmjNra2rgO/9/VzuFgWltbJem2vw75kdAIuFwuFRUVKSMjQ5mZmaqqqlJvb6+Ki4vD3dqo8Fd/9Vf62te+ppkzZ6qjo0OVlZWKiorSqlWrwt2a0S5cuBD0f1gffvihWltbNWXKFN19993asGGDvvvd72ru3LmaNWuWtmzZosTEROXn54evacMMdw6nTJmi73znO1q5cqXsdrvOnj2rb3/725ozZ45yc3PD2LU51q1bp9raWv385z/X5MmTA/elxMXFacKECYqLi9PatWvlcrk0ZcoUxcbG6umnn1Z2dra+9KUvhbl7M1ztHJ49e1a1tbX6yle+oqlTp+qdd97Rs88+q/vvv1+LFy8Oc/dhFu7HlEarF1980br77rut6OhoKzMz0zpx4kS4Wxo1CgsLrenTp1vR0dFWUlKSVVhYaP3qV78Kd1vGO3r0qCXpiq2oqMiyrM8ebd6yZYuVkJBg2Ww268EHH7TOnDkT3qYNM9w5vHjxovXQQw9Zf/RHf2SNHz/emjlzplVSUmJ5vd5wt22Mwc6dJOvll18O1Pz+97+3/uIv/sK66667rIkTJ1qPPvqo9fHHH4evacNc7Ry2tbVZ999/vzVlyhTLZrNZc+bMscrKyqzu7u7wNm6ACMuyrFsZkAAAAELFPSwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGO//AFPAfbl8GNdjAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(range(len(four_spec_8)), four_spec_8)\n",
    "plt.bar(range(len(four_spec_12)), four_spec_12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5.8845586355845466e-05+0j)"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(four_spec_8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
